{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 105,
   "id": "2a0b97d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.interpolate import interp2d, interpn\n",
    "from numpy.matlib import repmat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00f6be88",
   "metadata": {},
   "source": [
    "# Entry/Exit with Homogeneous firms\n",
    "## and CES demand\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "30a07c92",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sequence_jacobian import solved, simple, create_model\n",
    "\n",
    "@simple\n",
    "def household(w, psi, nu):\n",
    "    L = (w/psi)**nu\n",
    "    return L\n",
    "\n",
    "@solved(unknowns={'v': 1.0, 'N': 1.0}, targets=['valfunc', 'accumulation'])\n",
    "def firms(C, sigma, v, beta, delta, Ne, N, L):\n",
    "    mu      = sigma/(sigma-1)\n",
    "#    y       = p**(-sigma)*C\n",
    "    ell     = L/N            # imposing labor market clearing\n",
    "    pi      = (1-1/mu)*(C/N)\n",
    "    valfunc = v - (pi + beta*(1-delta)*v(1))\n",
    "    accumulation =  N - (1-delta)*N(-1) - Ne\n",
    "    return mu, ell, pi, valfunc, accumulation\n",
    "    \n",
    "@simple \n",
    "def mkt_clearing(mu, N, sigma, v, fE, w, Z, C, L, pi, Ne):\n",
    "    p       = mu*w/Z\n",
    "    rho     = N**(1/(sigma-1))\n",
    "    goods_mkt = p - rho                   # as in BGM\n",
    "    free_entry = v - fE*w/Z\n",
    "    budget     = C - w*L - N*pi\n",
    "    return goods_mkt, free_entry, budget\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e0a70a25",
   "metadata": {},
   "source": [
    "### Steady State"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "d3fbc144",
   "metadata": {},
   "outputs": [],
   "source": [
    "ces = create_model([household, mkt_clearing, firms], name=\"CES Model\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "bb885ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "calibration = {'beta':    0.96,\n",
    "               'psi':     1,\n",
    "               'fE':      1.0,\n",
    "                'nu':     2.0,\n",
    "                'delta':  0.11,\n",
    "                'sigma':  5,\n",
    "                'w':      1,\n",
    "                'Z':      1.0,\n",
    "                'Ne':     1.0, \n",
    "                'C':      6.5}\n",
    "\n",
    "ss = ces.steady_state(calibration)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "283c6fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'w':1.0, 'Ne': .25}, \n",
    "      targets={'free_entry': 0.0, 'goods_mkt': 0.0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "87cc0f0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ss=ces.solve_steady_state(calibration, \n",
    "      unknowns={'psi':1.0, 'Ne': .17, 'C': 1.07},\n",
    "      targets={'free_entry': 0.0, 'goods_mkt': 0.0, 'budget': 0.0})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e52b213",
   "metadata": {},
   "source": [
    "# Shock to Z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "id": "4595522d",
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 300\n",
    "dZ = .01 * .979 ** np.arange(T)\n",
    "irf = {}\n",
    "irf = ces.solve_impulse_linear(ss, ['w', 'Ne', 'C'], ['budget', 'free_entry', 'goods_mkt'], {'Z': dZ})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "id": "c214c308",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe5efdec880>]"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAD4CAYAAADo30HgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAX+klEQVR4nO3de5Bc1X3g8e+ve3pektCboOhhgVFggfViUAjZtUlsU0b2piKcNVuKqwyVOKvYZdcu5X0Eiq0tZ2upXSfrUOV4zUYJLsB5AGUWw8ambIjtkLhsQFqEHgjh4WWNNCChx+gxI83r7B99u6dbmhkBo57unv5+qrrmzrl9+57DRf2b33ncGyklJEkqydW7ApKkxmJgkCRVMTBIkqoYGCRJVQwMkqQqbfWuwHQtWbIkrV69ut7VkKSmsmXLlrdSSksn2tf0gWH16tVs3ry53tWQpKYSEa9Pts+uJElSFQODJKmKgUGSVMXAIEmqYmCQJFUxMEiSqhgYJElVDAyS1CCGRsZ48NmfMzpW38chNP0CN0maLe568iXu/tHLzO8qsO6KZXWrhxmDJDWIZ149BEBnIV/XehgYJKlB9Ow/DkBE1LUeBgZJahD9g8MADI+M1bUeBgZJagCHTgyVt4dHDQyS1PJefONoeXvIwCBJOjIwXN4eHq3vdFUDgyQ1gIGh0fK2XUmSJAaHRsrbBgZJ0mkZg11JktTyTtiVJEmqNDg0Qnu++JXsOgZJEgNDo8zrbCPCjEGSBAwOjdLVnqeQzzHkGIMk6cTQCHPa22jP58wYJGm223NogG/846ukNHkmMFDOGMLAIEmz3be29PJf//YF9hwarCofHUvs2NsPFLuSurOuJAODJM1yvYeLAWHb3iP0Dw7z7ef2klLi+zvf4Df+9B/5+cEBBioCw9CIYwySNKvc++NX+bU//iGHszum7j0yAMD23n6+/dxebn1wKy8fOM6rB08AsK9/kIGhEbrb2+xKmkhErIuI3RHRExG31bs+kvROpJT4y6d/zusHB/jS/90JwN4jWcbQ209f/0kAXug7xpvZ9qETQ1UZg4GhQkTkgf8FfAy4DPjtiLisvrWS1MoGhkYYqLiPERS//Pv6Bzl4/FRV+dhY4tnXDtOz/zhrzp/Lo1v3seX1Q/QdKQaAHXv7eaO/GCR29R0tB4mDJ4aqpqvWOzC01fXsZ7oG6EkpvQIQEQ8A64EXzvWJ9hwa4JW3Tpzrj31bppqZMCPnb9mTQ6pjBep52ev8v1xdL/tYShwZGGL14jnM6WjjyMAwRwaHKORzzOto42f7j5PLBVevWkhbPtje208ELOguML+rwL9/6HmODA7zh795OZ2FPEMjY+zY28+fPfUKS+d18L1br2Nhd4Gh0TG+8v2X2PTUK+QC7rnll7nx6z/mtoe3MzKW+KfL57N9bz/P7TkCwIt9R8sP5zl0fIiB4SxjaKv/OoZGCwzLgT0Vv/cCv3L6myJiI7ARYNWqVe/qRN/d3sd/f/zFd3WspNbR0Zbj0mXnceuDWwmgkM8RAddetIgtrx/mtoe30d2e59nXDnN0cJhrL1rE733gIlYt7mbjdRfxP7LvmQ+sWcL2vf28frA43rCr7xhjWcR+4+ggo2OJ7vY22vPBiBlDlYmegH1G6EwpbQI2Aaxdu/Zdhdb1Vy5n7epF7+bQc6LOz/qe8D/0jJ27zo2vb9vreO66try+bT+vs8DuN4tfxAu725nfVWB4dIyjg8OsXNRNLhf8uOct8hG8b8V8IuDYyRFeevMYqxfP4Z+tXMDvf3MLhweG2NZbnF76B+suZcvrh/lv39lVda7/eMOlXP2ehQB86ldWjQeGi5dw949eBqCQD944erJ8TGkaa6OMMTRaYOgFVlb8vgLYV4sTXTC/kwvmd9bioyU1oFWLu6fc/6/Xrjyj7P2rFpa37/vda0gp8W/u30xf/0muXLmAK1cu4OTwKD/bf5zr1izl+d4jXLVqQfmY8zoLfOiSpfxw9wHev2oBnYUcJ4fH+OjlF/CdbX3l9+05XMwiSoGh8k6r9dBogeFZYE1EXAjsBTYAn6pvlSSpKCL42qeuYiylcub7hQ+vKe//V1evOOOYTTevZe/hQbrb21i9eA4vvnGMD19yfjkwdBZy5XUOXe1txYzBu6uOSymNAF8AvgfsAh5KKe2sb60kaVxnIU93+9v/m7qQz7F6yRwAVi8u/rxgfifvzzKLi8+fy+hYsUd8Tnue9rb6r2NotIyBlNJ3ge/Wux6SdK69Z0mxO+sXzuvgzz59NX/99M85OTzGjr1HARpmumpDZQySNJt94OIlrDl/LisWdnP+vE5uvf6XWDqvo7y/u9SVlE1Xff3gCf5u15szPkup4TIGSZqtPrhmKU988deqypbMbQfg/HkdXHrBPAr5YGh0jJ79x7j+T54C4N7f+WV+/ZLzZ6yeBgZJqqNfv+R8/vO//Cf89jWr6CyMdyXtPzq+qvrwwNCM1snAIEl1NL+rwO998KLy76VZSSdHxqesHjs5MtGhNeMYgyQ1kNIYw+DQ+LiCgUGSWlh7NsYwOGzGIEmimDEAnDg1HgyOnxqe0ToYGCSpgRTail/Lx04Wg8GiOe0cN2OQpNZVyhiOZsFgydx2u5IkqZW154v3YDp2cph8LljY3c6xUwYGSWpZ5YxhcISuQp55nQWOnRzhp68cnLEV0AYGSWog411Jw3QWcszrbGNX31E2bPopm/7hlRmpg4FBkhpIafD56MkROgt55naMr0N+o//kZIedUwYGSWoghVw2xjA4TGchz7zO8cCweE7HZIedUwYGSWoglbOSugp55lYEhva2mfnKNjBIUgMZ70rKxhgqupJODs/MIz8NDJLUQDqzwDA0MpZ1JRXK+wYNDJLUehbNaS9vd502+Dw4ZGCQpJZTGRg6C3nGUir/bsYgSS1oQXc7UZyYRFchzy8u6CrvMzBIUgvK54IFXcVxhc5CjiuWz+cf/tOHuPwXz+OkXUmS1JpK3Umd7XkAVi7qprs9z4CBQZJaUzkwtOXLZZ2FvF1JktSqutuLM5G62vMVZXnXMUhSq+rOAkJXYTwwdJkxSFLrKmUKlbfA6HKMQZJa15ysK+lURYbQWcg7K0mSWlWpK2mgIjB0t9uVJEkta+m84u212/MVXUmFPCNjieEZeIpb29nfIkmaSbf889WMpcSnf/U95bLObCB6YGiU+V21/ZvejEGSGkwhn2Pjde+lo61yumrx7/iZmLJas8AQEV+KiL0RsTV7fbxi3+0R0RMRuyPihoryqyNie7bvqxGlO4ZIUmvrai9+Xc/EHVZrnTHclVK6Mnt9FyAiLgM2AJcD64CvR0QpLN4NbATWZK91Na6fJDWF0pqGmRiArkdX0nrggZTSqZTSq0APcE1ELAPOSyn9JKWUgPuBG+tQP0lqOJVjDLVW68DwhYjYFhHfiIiFWdlyYE/Fe3qzsuXZ9unlZ4iIjRGxOSI2HzhwoBb1lqSG0jRjDBHxZETsmOC1nmK30HuBK4E+4Culwyb4qDRF+ZmFKW1KKa1NKa1dunTpdJogSU2h3JU0AxnDtKarppSufzvvi4g/B/42+7UXWFmxewWwLytfMUG5JLW80uDziaGRmp+rlrOSllX8+glgR7b9GLAhIjoi4kKKg8zPpJT6gGMRcW02G+lm4NFa1U+SmsmC7uKtuPsHh2t+rloucPujiLiSYnfQa8DvA6SUdkbEQ8ALwAjw+ZRSKTf6HHAv0AU8nr0kqeUt7G4nF/DWsVM1P1fNAkNK6dNT7LsTuHOC8s3AFbWqkyQ1q3wuWDSnnbdODNX8XK58lqQmsXhOx4xkDAYGSWoSi+e2c9CMQZJUsmRuBwePn+LEqRH+5/d2s3XPkZqcx8AgSU1i8dx23jo+xIlTI3zthz3s3Ndfk/MYGCSpSSyZ28HxUyPl22LkanSfUQODJDWJJXOLaxneOl4cgM7V6P7TBgZJahKL5xSf7LY/m5lUqycTGBgkqUkszjKGA8dKGYOBQZJa2qI5xcBQmrJqV5Iktbh8FglGRscAMwZJanmlQDA6VnwiQa0efmxgkKQmUQoMI1lgMGOQpBZXGlMYNTBIkmB8eupweYyhNucxMEhSkzg9Y3AdgyS1uPKspHJXUm3OY2CQpCYRp81KcoxBklpcKUMoZww1+gY3MEhSkxhfx1AcfHaMQZJaXK48K8muJEkS4yudR5yuKkmCM1c+580YJKm1laaruo5BkgRMMCvJriRJam1nrGOoUWQwMEhSE8mFg8+SpAq5iHJXkmMMkiRyEd4SQ5I0Lpdz8FmSVMGMQZJUpTjGULpXUo3OMZ2DI+KmiNgZEWMRsfa0fbdHRE9E7I6IGyrKr46I7dm+r0Y2ehIRHRHxYFb+dESsnk7dJGk2ioDRBr9X0g7gt4CnKgsj4jJgA3A5sA74ekTks913AxuBNdlrXVb+GeBwSuli4C7gy9OsmyTNOrkIhhu5KymltCultHuCXeuBB1JKp1JKrwI9wDURsQw4L6X0k5RSAu4Hbqw45r5s+1vAR6JWc7EkqUnlovJBPTU6R20+luXAnorfe7Oy5dn26eVVx6SURoB+YPFEHx4RGyNic0RsPnDgwDmuuiQ1rlxEeYFb3dYxRMSTEbFjgtf6qQ6boCxNUT7VMWcWprQppbQ2pbR26dKlUzdAkmaRXC5qnjG0ne0NKaXr38Xn9gIrK35fAezLyldMUF55TG9EtAHzgUPv4tySNGvlonIdQwOOMUzhMWBDNtPoQoqDzM+klPqAYxFxbTZ+cDPwaMUxt2TbnwR+kI1DSJIyM7GO4awZw1Qi4hPAnwJLge9ExNaU0g0ppZ0R8RDwAjACfD6lNJod9jngXqALeDx7AdwDfDMieihmChumUzdJmo0q75WUq9Gf9tMKDCmlR4BHJtl3J3DnBOWbgSsmKD8J3DSd+kjSbFeZJDRbV5IkqQYqg4GBQZJUNROp2dYxSJJqoPKpbT6PQZJ0WldSjc5Rm4+VJNVCzsFnSVIlB58lSVUqxxWiRt/gBgZJaiJ2JUmSqjj4LEmqUjld1YxBklSVJTTkM58lSTPLWUmSpCoOPkuSqoSDz5KkStVjDGYMktTy8llkqFW2AAYGSWoqpXGFfA0jg4FBkppIqfuoVt1IYGCQpKZSShTsSpIkAeNdSbWaqgoGBklqKuMZg4FBkkTlGEPtzmFgkKQmkrcrSZJUKZd9azv4LEkCxruSzBgkScB4QHAdgyQJcB2DJOk0rmOQJFUJMwZJUqW8YwySpErlrqQafntP66Mj4qaI2BkRYxGxtqJ8dUQMRsTW7PW/K/ZdHRHbI6InIr4aWdiLiI6IeDArfzoiVk+nbpI0G42vY2jcjGEH8FvAUxPsezmldGX2+mxF+d3ARmBN9lqXlX8GOJxSuhi4C/jyNOsmSbNOw69jSCntSintfrvvj4hlwHkppZ+klBJwP3Bjtns9cF+2/S3gI1HLTjRJakKlQedmvVfShRHxXET8fUR8MCtbDvRWvKc3Kyvt2wOQUhoB+oHFE31wRGyMiM0RsfnAgQO1qb0kNaDyE9xqGBnazvaGiHgSuGCCXXeklB6d5LA+YFVK6WBEXA18OyIuByZqSSqdaop91YUpbQI2Aaxdu3bC90jSbDQT6xjOGhhSSte/0w9NKZ0CTmXbWyLiZeCXKGYIKyreugLYl233AiuB3ohoA+YDh97puSVpNss16223I2JpROSz7YsoDjK/klLqA45FxLXZ+MHNQCnreAy4Jdv+JPCDbBxCkpRp+Af1RMQnIqIX+FXgOxHxvWzXdcC2iHie4kDyZ1NKpb/+Pwf8BdADvAw8npXfAyyOiB7gi8Bt06mbJM1GuVzt1zGctStpKimlR4BHJih/GHh4kmM2A1dMUH4SuGk69ZGk2S4aPWOQJM0sb7stSaribbclSVV85rMkqcr4LTFqdw4DgyQ1EccYJElVHGOQJFUpr2MwY5AkgesYJEmnadp7JUmSasPpqpKkKuHgsySp0kw8j8HAIElNpDxdtYYpg4FBkprI+HTVGp6jdh8tSTrXwq4kSVKlhn+CmyRpZuVdxyBJquSsJElSFdcxSJKqmDFIkqrksm9tn8cgSQIqM4YanqN2Hy1JOtfsSpIkVSkHhhp+exsYJKmJlLqQHGOQJAGVt8So3TkMDJLURLwlhiSpioPPkqQq4+sYaniO2n20JOlcM2OQJFUpBYR8oz7BLSL+OCJejIhtEfFIRCyo2Hd7RPRExO6IuKGi/OqI2J7t+2pkQ+wR0RERD2blT0fE6unUTZJmo1wT3Hb7CeCKlNL7gJeA2wEi4jJgA3A5sA74ekTks2PuBjYCa7LXuqz8M8DhlNLFwF3Al6dZN0madRp+VlJK6fsppZHs158CK7Lt9cADKaVTKaVXgR7gmohYBpyXUvpJSikB9wM3VhxzX7b9LeAjUcsVHJLUhJptHcPvAo9n28uBPRX7erOy5dn26eVVx2TBph9YPNGJImJjRGyOiM0HDhw4Zw2QpEY3ExlD29neEBFPAhdMsOuOlNKj2XvuAEaAvyodNsH70xTlUx1zZmFKm4BNAGvXrp3wPZI0G+VypTGGOgaGlNL1U+2PiFuA3wA+knUPQTETWFnxthXAvqx8xQTllcf0RkQbMB849DbaIEkto+Fvux0R64A/AH4zpTRQsesxYEM20+hCioPMz6SU+oBjEXFtNn5wM/BoxTG3ZNufBH5QEWgkSTRIV9JZfA3oAJ7I0pqfppQ+m1LaGREPAS9Q7GL6fEppNDvmc8C9QBfFMYnSuMQ9wDcjoodiprBhmnWTpFlnJjKGaQWGbGrpZPvuBO6coHwzcMUE5SeBm6ZTH0ma7cbXMTTodFVJ0syKRl/HIEmaWQ0/+CxJmlmlu6uaMUiSAMg3wb2SJEkzKLzttiSp0vg6hhqeo3YfLUk618qDz436PAZJ0sxyHYMkqUopHuQNDJIkGH+kp2MMkiQAls3v5N9++GI+dOn5NTvHdG+iJ0maQRHBFz96SU3PYcYgSapiYJAkVTEwSJKqGBgkSVUMDJKkKgYGSVIVA4MkqYqBQZJUJVJK9a7DtETEAeD1d3n4EuCtc1iderItjcm2NCbbAu9JKS2daEfTB4bpiIjNKaW19a7HuWBbGpNtaUy2ZWp2JUmSqhgYJElVWj0wbKp3Bc4h29KYbEtjsi1TaOkxBknSmVo9Y5AkncbAIEmq0rKBISLWRcTuiOiJiNvqXZ93KiJei4jtEbE1IjZnZYsi4omI+Fn2c2G96zmRiPhGROyPiB0VZZPWPSJuz67T7oi4oT61ntgkbflSROzNrs3WiPh4xb6GbEtErIyIH0bErojYGRH/LitvuusyRVua8bp0RsQzEfF81pY/zMpre11SSi33AvLAy8BFQDvwPHBZvev1DtvwGrDktLI/Am7Ltm8Dvlzvek5S9+uAq4AdZ6s7cFl2fTqAC7Prlq93G87Sli8B/2GC9zZsW4BlwFXZ9jzgpay+TXddpmhLM16XAOZm2wXgaeDaWl+XVs0YrgF6UkqvpJSGgAeA9XWu07mwHrgv274PuLF+VZlcSukp4NBpxZPVfT3wQErpVErpVaCH4vVrCJO0ZTIN25aUUl9K6f9l28eAXcBymvC6TNGWyTRyW1JK6Xj2ayF7JWp8XVo1MCwH9lT83svU/+M0ogR8PyK2RMTGrOwXUkp9UPzHAdTuaeHn3mR1b9Zr9YWI2JZ1NZXS/KZoS0SsBt5P8a/Tpr4up7UFmvC6REQ+IrYC+4EnUko1vy6tGhhigrJmm7f7L1JKVwEfAz4fEdfVu0I10ozX6m7gvcCVQB/wlay84dsSEXOBh4FbU0pHp3rrBGWN3pamvC4ppdGU0pXACuCaiLhiirefk7a0amDoBVZW/L4C2FenurwrKaV92c/9wCMU08U3I2IZQPZzf/1q+I5NVvemu1YppTezf8xjwJ8znso3dFsiokDxi/SvUkr/JytuyusyUVua9bqUpJSOAD8C1lHj69KqgeFZYE1EXBgR7cAG4LE61+lti4g5ETGvtA18FNhBsQ23ZG+7BXi0PjV8Vyar+2PAhojoiIgLgTXAM3Wo39tW+geb+QTFawMN3JaICOAeYFdK6U8qdjXddZmsLU16XZZGxIJsuwu4HniRWl+Xeo+61+sFfJzibIWXgTvqXZ93WPeLKM48eB7YWao/sBj4O+Bn2c9F9a7rJPX/G4qp/DDFv3A+M1XdgTuy67Qb+Fi96/822vJNYDuwLfuHuqzR2wJ8gGKXwzZga/b6eDNelyna0ozX5X3Ac1mddwD/JSuv6XXxlhiSpCqt2pUkSZqEgUGSVMXAIEmqYmCQJFUxMEiSqhgYJElVDAySpCr/H+zLAI9BI7w5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(irf['L'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "211acf94",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-113-287d6a2f06e0>:4: RuntimeWarning: invalid value encountered in power\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe5efea7fd0>]"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAjFElEQVR4nO3dfZRcdZ3n8fe3qro76aTJc0LniRCMaBAJ2AMqiiKEgThjwBlm0Vkmu0eNcxbOqDseNx7mOLh7jssyPozOcWWDonHUYRw1EiXKQw4eCMpDI0lICCEhBPLQJE0eSNJ56If67h91bz11VT/UrVB9uz6vc+rUvb/7u1W/Wze53/493WvujoiI1K9ErQsgIiK1pUAgIlLnFAhEROqcAoGISJ1TIBARqXOpWhegElOnTvV58+bVuhgiIrHyzDPPvO7u04rTYxkI5s2bR3t7e62LISISK2b2Sql0NQ2JiNQ5BQIRkTqnQCAiUucUCERE6pwCgYhInVMgEBGpcwoEIiJ1ToFARGKrty/NT9t305fW7fSjqEogMLNrzWybme0wsxUltv+1mW0KXr83s4vytu0ys+fMbIOZaZaYiAxZ+yuH+cLPNvHsq4drXZRYizyz2MySwLeBxcAe4GkzW+Puz+dlexn4gLsfNrPrgJXAZXnbr3T316OWRUTqS09fOnhXjSCKatQILgV2uPtOd+8G7gWW5mdw99+7exiynwBmV+F7RaTOhQ9Y1JMWo6lGIJgF7M5b3xOklfMJ4Dd56w48aGbPmNnycjuZ2XIzazez9s7OzkgFFpHRIbz8F3cRvHKwS8FhGKoRCKxEWskzYGZXkgkE/yMv+XJ3vwS4DrjFzK4ota+7r3T3Nndvmzat383zRKQOhRf7dN5Ff++Rk3zwq79j/Q61Ng9VNQLBHmBO3vpsYF9xJjN7J/BdYKm7HwzT3X1f8H4AWE2mqUlEZFDh9T8/EBw50Y07HOrqrlGp4qcageBpYIGZnWtmjcBNwJr8DGY2F/gFcLO7v5iXPs7MWsJl4BpgcxXKJCJ1wIPGh/xWoHAoqTqQhy7yqCF37zWzW4EHgCRwj7tvMbO/DbbfBXwJmAL8XzMD6HX3NmAGsDpISwE/cfffRi2TiNSHUjWCMAD0BiOKZHBVeTCNu68F1hal3ZW3/EngkyX22wlcVJwuIjIUaS98h7wagSaZDZlmFotIbJXqLO5NZ2oCqhEMnQKBiMRWePnPHyram20aUo1gqBQIRCS2cjWCXFquaUg1gqFSIBCR2CrVWdybVo1guBQIRCS2Ss0sDvsGetVZPGQKBCISW2FNwEvWCNQ0NFQKBCISW6WahsI+AtUIhk6BQERiK9s0lPfHf+7W1IU1glM9fcxbcT8/euKVYX3HqZ4+untHd+1CgUBEYqvUPIK+Mp3FR0/2APDPD7/IcPyX7z/Fl3+1JUoxR7yqzCwWEamF3PMIcmnZPoKi4aNhlhPdfcP6jp2dXYz2O1qrRiAisRXedK5g+GiZp5aFTUX5gWDz3jc4eqpnwO84crKHg6P8TqYKBCISW+Ef/elSNYKiPoLipiJ358/+ZT3/+btPFqT/z189z5qNmTvph/0DB4+frnLJRxY1DYlIbOXmEfTvIyi+6VxxU1GYb9OeNwrS73n8ZXgcPnLRTI6cyNQWjpzsobcvTSo5Ov92Hp1HJSJ1wYcxj6B/U1H/hv9TPYX9B0dOdgefD4dPDNyEFGcKBCISW7l5BLm0cjedK14vdS+i/L6AQ13d2RpBZtvobR5S05CIxFapzuK+4AJf3DRUfOHvKTE34NDxXCB44bWjHDvVm10/eHz0dhgrEIhIbJWqEfSUaxrKu/DvPnSi5OcdOpEXCDqOMb4pd4l8fZAO48Nd3Zzo6WPWxLFDKvtIUpWmITO71sy2mdkOM1tRYruZ2beC7ZvM7JKh7isiUk4YAHwIE8rybznx/jsf4bZf5h6Png62Hcpr/vmnB7bxg9/vyq4fPN7Nmo37WLd1f79ybN9/jGu/+Sh//i/rOX66t9/2kS5yIDCzJPBt4DpgIfAxM1tYlO06YEHwWg58Zxj7ioiUVHoeQekJZcW3nHhy58Hs8hvBrOOw+edfP3EpZ08Yw/MdR0kljGTC+O2W1/jsvc/ypfu2ZAPPye4+3jjRw60/eZbTvWkOdXVz96M7CwLTcPT2pXnlYFfF+1eqGk1DlwI7gucPY2b3AkuB5/PyLAV+6Jmje8LMJppZKzBvCPuKiJRUsrM4Xfo21MU1hNYJY9h1MNNEdOhEN5PGNXKoq5tUwnjfW6by/gVTefn1Ls4a28Chrm6eevkQDUlj75GTPLb9dczgH+/bwp4jJ+nuTfPtj1/C6mf38s1123lgy2t89JJZtE4YS9u8SXxr3XZ2HDhOMmEcOHqa5qYkLU0NpJJG57HTnOzpY/aksew+dJJXD51g/tRxTGhuyByXO2mHRMK4cNZZLHvPPBbMaKnq71iNQDAL2J23vge4bAh5Zg1xXwDMbDmZ2gRz586NVmIRGRVKP7M4mEfQr2mosEYwraUpGwgOd3XDNDgcBAQz44KZZwGZ2sK8Kc3sOniCn3zq3dx41x/4m3ueAqClKcVFsyfQ3JhiyYVnc+XbprFmwz6+/bsdfGXtCyQTxrvmTmLjniO8c/YETvWkefvMszjZ3cexUz2c7HFmT2qmqSHB3sMnOWdKM//pT+bw9K5D9KWdhBlmkDDjZHcfq/+4lxsunlX137EagcBKpBXXa8rlGcq+mUT3lcBKgLa2tlF+5w8RGYrcM4tzaX3Z4aNDn0fw2tFTrH52D8++eoQp4xoBuGDmhMznpZ0ff+rdNKUSTB3fxJXnT2PzvqN85YYLecess2idkOscbm5McdOlc/mLd81m4+4j/OVdf+CpXYf41PvP5bYPR2/17ulLk7BSl81oqhEI9gBz8tZnA/uGmKdxCPuKiJSUbRrKawbqKdc0VFQjONGd69S99SfPZpffM38KAAtmjM+m5Y8Euuvmd5FKJEgmyl+QG5IJ2uZNZtGciWzYfYQb2+aUzTscDWdoZnM1PvVpYIGZnWtmjcBNwJqiPGuAvwlGD70beMPdO4a4r4hISels01AuLXuLiUFqBF2nC2cR//lFMwF4NRha2pRKsmjORD59xfyCfE2p5IBBIN8Xrj2fz1y1gLdWuU2/2iLXCNy918xuBR4AksA97r7FzP422H4XsBZYAuwATgD/daB9o5ZJROrDcB5eH65/8Pxp/G5bZ0GN4GOXzuUfPvx2frVxH1e8dWo2/Ze3XB6pfO89byrvPW/q4BlrrCoTytx9LZmLfX7aXXnLDtwy1H1FRIYi10fQ/zbU+U1BW/a9kb3wf/XGi/j43U9kO4p/ecvlLJozEYCNX7qGsY3JM1/wEUYzi0UktnzApqHM+7bXjvHhb61nXHCBb0gkaEgmso+fbMxrd5/Q3PBmFHvE0U3nRCS2Bm4aylzon3w5M3GsK3ggTSppBZ2uDcnqj8KJGwUCEYmt3MziXFrYFxDec+j5fUcBGNsQ1AiSCRpT+YFAl0H9AiISW7lnFufXCII+gqBG8HxHJhCc6s3UCBqSVtAclFKNQIFAROIrXaJpKOwjSHvmXkBbg0DgDsmEYWYFzUGNqhEoEIhIfJVqGsqfL/Di/mMF66lg/H9+09BoffzkcGjUkIjEVn5ncU9fmh898UrBbaB3HDhekD/sD1BncSEFAhGJrdwzi+Hux3Zy52+3FWx/qTMTCFqaUhw73ZvtD2hMqrM4n34BEYmt/BrBIy8c6Lf9pc7jTGpuYFpLEwCpROaSp1FDhfQLiEhsha3/fWnn6V2H+21/qbOLc6aMy174w2ag8OJvxpDvGzSaKRCISGyFo4Ve3H+s5PYdB44zb0ozTcEcglRRIFBtIEO/gojEVtg0dPhET0F6Ku+v/LlTxtEU1giKmoYaVBsAFAhEJMbCpqHiW06PacjdOG7u5OZsIMh1Fgc1g5QugaBAICIxFo4aKr7l9JiG3KWtdcIYmlJB01CisEkoXK93+hVEJDbcnd3Bg2My65n34hpBeOEHmHHWGJqCwBDWAML3Rs0hABQIRCRG1u94nSv+6RH2HjkJ5GYWFz+WMr9GcPaEMYxJhbegLpxHoFnFGfoVRCQ2Oo+dxh0Od3UDuVtLFNcI8h8uM74pla0RpIr6BjSrOEOBQERiI3yYzOngTqJh01DYR9AyJnOzhDGpwqeMZUcNJQubhDR8NCPSr2Bmk83sITPbHrxPKpFnjpk9YmZbzWyLmX0mb9vtZrbXzDYEryVRyiMio9vpbCDIvIdNQ2GNoKUpCAQNxYEg7CwuvOmcAkFG1F9hBbDO3RcA64L1Yr3A37v724F3A7eY2cK87d9w90XBS88uFpGywppANhCENYKgjShsEgoDwYSxmUdP5oaPFo0aUtMQED0QLAVWBcurgOuLM7h7h7v/MVg+BmwFZkX8XhGpQ9mmoZ4wEBR2EoeBILzwTxnfmFlvKOwT0MziQlF/hRnu3gGZCz4wfaDMZjYPuBh4Mi/5VjPbZGb3lGpaytt3uZm1m1l7Z2dnxGKLSJy8caKH7z/+MqeCANDdV1gjCIWPozzRnbkV9dRxmZvNFc8jKL73UL0bNBCY2cNmtrnEa+lwvsjMxgM/Bz7r7keD5O8A5wGLgA7ga+X2d/eV7t7m7m3Tpk0bzleLSMx95t+f5cu/ep6Ne44AcLon00SULg4EjZk+ggtmTgDg765aAFBiZrFqBPkGfR6Bu19dbpuZ7TezVnfvMLNWoP99YDP5GsgEgR+7+y/yPnt/Xp67gV8Pp/AiUh+effUIkOsbyNYIKIwEzUGNYGJzA7vu+HA2vfheQ2oaKhT1V1gDLAuWlwH3FWcwMwO+B2x1968XbWvNW70B2ByxPCIyCr1xMnNTuaPBe66PoDBf2EcQBoxQ8d1H1TRUKGoguANYbGbbgcXBOmY208zCEUCXAzcDHyoxTPROM3vOzDYBVwKfi1geERll0nntP9mAcKqH7z62M/ug+lA4WqhfIEgV1wQ0jyBfpEdVuvtB4KoS6fuAJcHyeqBk2HX3m6N8/3D19KVJmOlBFCI19v3HX+bep3ZzwyWz+PQV88k0HJS25/DJ7HIYCP754e0ATGpuKMgbdhZ3lw0ERbeY0E3ngDp7ZvH/+vXz/PAPr5BMGI3JBI2p4JVM0JRK0FCUll1OJWhKFm3Py9OQNJKJBKmEkUpa5j2RIJXMBJ1UwbZENk8yYTQkE8F74WckE0ZDIkEyyJuw8JV5otJA/3FERopnXz3MsVO9XPHWaRw50c3Jnj6aG1J8/cEXSSWNO37zAj97Zg/jm1Ks/m/vxcxIp52u7l5axjSwee8bfOFnm7Kfd6K7r+Dze/qKh49mLuz9A0HYNFQ4aqgxpf9HUGeB4Mq3TWfyuEZ6+tJ09wavvjSne3PrPX2ZtO7eNCdO9Ga29RVtD9KK/xG+2ZIJI2mWfdxeGCgSQXoiEQQOywSOZCK3PWFFeYLgkjQyn5P9DDAy32FmGJnH+2XeM/sSbs9uy+1HUd5wOdgtk9dy+1m4XPBdVvCdxWmJROHnhekE+4csr2IaphdcBoJE659UsH9hWol8JYL0UD+nsGz9y1v4mXnbS35PiXwltqc996SvgvNqxrjGJC1jGjh4/DQTmxsY15TiZ8/sAeC6d5zN2MYUHUdOsnDmWfy0fTe/f+kgsyaO5eI5E+lNO9997GW6+9J8+MJW1u94nYTBdRe2cux0L/f/3fv4weO7+I/g8w51dTNlfBN3P7aT//2bF3jmH67mR0+8wkudx3nb2S288Fr/p5A1phJwOreerRH0FQaM7DyCRGGTkGoEGfUVCM6fzpXnDzjVYVjSaae7L01f2ulNO73Bck/a6etzetLBtj6nN50O8gTLfV6wX2+6MD3zGZn0vrTT54575tmsaXfSQVraM+VIu9OXzvyHziwXbcvbP7NcKg/Z9PDV3ZfJ62Q65gqWyXxm2nMTewrTgzEdTrBP5vvC7WFHX7/0ft9V9P15nxemF35Xriyh/JCdLWtBWtX+WYx6LU0pGlMJ7tuwr9+2D71tOq8c7GL9jtdxh8vOnczcyc387sVO3t7awlMvH+InT77KX7XN5oKZE7jzL9/JNReczad+2M4Pfr+LrR3HeP145sp+/3MdPLx1P1cvnME//tlCLv3Kun7f11T0YJkp4zPzBpobUyXzFc8sVh9BRl0FgmpLJIwxieTgGSWW8metZoNWie3lAkqpgDRQvkq/h+w+Q/uc/HR3ghqZZfPkB92u7j6OnuxhUnMjh090c/x0LxfMPIvxTSkef+kgBsybMo6ftu+mdeIY/vqyc4DcHyzFF9qvP7iNP+w8yO0fuQDIfO9508YBcPdjOznVk2Z+sP7Nh7dzsKubaxbOyI76KdZYFAiWLppJ1+lePn7Z3IL0XNNQ8TwCNQ2BAoFIWeWaVPJyvGllGQnmTmkuWP/AW3MTOz//p+cXbEsmjGSJ3+e/X3N+v7TZk5pJGNlZwzs7uwA42NXNrIlj+eD507NNPsWSRScmlUjwyffP75ev7DOLVSMAFAhEpMYaUwlmTRrL7kO50UGf/sB8bnzXbOZOHkdjKoG7k7D+M4n7iqpb5UYEhsNKc8NGMwM0mptUowcFAhEZAeZNGVcQCOZPHcdbprdk182MsQ1JuopGDRXPIyg3Mnzq+Ea+eN3buPYdmTmsqWSCH33iMs4/u6X0DnVG9SIRqbl5UzL9AmH/wJzJzf3yFD9jAAoDQX5fRzEz49MfOI+zJ4zJpl02fwoTmxsjlXu0UI1ARGrur9rmMGV8I53HTrOzs4tzgsCQb7BAoImilVMgEJGau3D2BC6cPYH1219n/9HTtJ41pl+e/AfSh9KeXyNQIKiUAoGIjBjvWzCV9y2YWnJb/gPpQ6oRVIf6CEQkFoofSA8KBNWiQCAisVCqRpA/aEiBoHIKBCISC02D1QjUR1AxBQIRiYWwRpB/W4mC4aOqEVRMgUBEYmFMEADOGpMb45I/s1g1gsopEIhILIQ1gpYxuYfRqLO4OiIFAjObbGYPmdn24H1SmXy7gkdSbjCz9uHuLyISTihrGVN61LseLVC5qD/dCmCduy8A1gXr5Vzp7ovcva3C/UWkjg0WCPSQmcpF/eWWAquC5VXA9W/y/iJSJ8KZxS1NDSW3q2WoclEDwQx37wAI3ss9/suBB83sGTNbXsH+mNlyM2s3s/bOzs6IxRaRuBk7SI1AfQSVG/QWE2b2MHB2iU23DeN7Lnf3fWY2HXjIzF5w90eHsT/uvhJYCdDW1qYHC4rUmVzTULkagQJBpQYNBO5+dbltZrbfzFrdvcPMWoEDZT5jX/B+wMxWA5cCjwJD2l9EZMLYTACYMr70raNVI6hc1KahNcCyYHkZcF9xBjMbZ2Yt4TJwDbB5qPuLiAAsXjiDH3/yMt573pSS2xUIKhc1ENwBLDaz7cDiYB0zm2lma4M8M4D1ZrYReAq4391/O9D+IiLFGpIJLn/LVC6eO4lHPv/BftvVNFS5SLehdveDwFUl0vcBS4LlncBFw9lfRGQgzSVuQKcaQeU08FZEYqfUH/8KBJVTIBCR2CnVDKR7DVVOgUBEYqdkIFCNoGIKBCISO6Wu+boNdeUUCEQkdqxk01ANCjJKKBCISOyU+uNfTUOVUyAQkdgp1UegeQSVUyAQkdhRZ3F1KRCISOxoHkF1KRCISOwoEFSXAoGIxI4mlFWXAoGIxE7JzmLVCCqmQCAisVNy+KhqBBVTIBCR2Ck1oUw1gsopEIhILBVf95O6mlVMP52IxFJxP4GahiqnQCAisdQvECR0OatUpF/OzCab2UNmtj14n1Qiz/lmtiHvddTMPhtsu93M9uZtWxKlPCJSP4orAGoaqlzUn24FsM7dFwDrgvUC7r7N3Re5+yLgXcAJYHVelm+E2919bfH+IiKlFNcI1FlcuaiBYCmwKlheBVw/SP6rgJfc/ZWI3ysida5fZ7H6CCoWNRDMcPcOgOB9+iD5bwL+rSjtVjPbZGb3lGpaEhEppXgIqW4xUblBA4GZPWxmm0u8lg7ni8ysEfgI8B95yd8BzgMWAR3A1wbYf7mZtZtZe2dn53C+WkRGoeIKgG5DXbnUYBnc/epy28xsv5m1unuHmbUCBwb4qOuAP7r7/rzPzi6b2d3Arwcox0pgJUBbW5sPVm4RGd36jxpSIKhU1KahNcCyYHkZcN8AeT9GUbNQEDxCNwCbI5ZHROpE/wllCgSVihoI7gAWm9l2YHGwjpnNNLPsCCAzaw62/6Jo/zvN7Dkz2wRcCXwuYnlEpE6oRlA9gzYNDcTdD5IZCVScvg9Ykrd+AphSIt/NUb5fROpXv85i9RFUTFMwRCSWiisAmkdQOQUCEYml/vcaqlFBRgEFAhGJpbACEPYNqI+gcgoEIhJLYR9B2DegpqHKKRCISCyFLUPhTUfVWVw5BQIRiaWwjyAVRALVCCqnQCAisRRe98P3lAJBxRQIRCSWwhqBOoujUyAQkViy7KihoGlIfQQVUyAQkVjK1QgI3hUIKqVAICKxlCgePqoaQcUUCEQklrJNQ0n1EUSlQCAisVRcI9DD6yunn05EYimcSJbIjhrS5axS+uVEJJaMcEJZYc1Ahi/S8whERGol7BKYO3kcR0/2cs6U5toWKMYUCEQklsKbzs2eNJbvLmurcWniTU1DIhJLYY1ALULRRQoEZnajmW0xs7SZlQ3JZnatmW0zsx1mtiIvfbKZPWRm24P3SVHKIyL1Ixw1FPYVSOWi1gg2Ax8FHi2XwcySwLeB64CFwMfMbGGweQWwzt0XAOuCdRGRQWUDgeJAZJECgbtvdfdtg2S7FNjh7jvdvRu4F1gabFsKrAqWVwHXRymPiNSP7PMIFAgiezP6CGYBu/PW9wRpADPcvQMgeJ9e7kPMbLmZtZtZe2dn5xkrrIjEQ65GoEgQ1aCjhszsYeDsEptuc/f7hvAdpc6SD2G/wh3cVwIrAdra2oa9v4iMLuH8MYWB6AYNBO5+dcTv2APMyVufDewLlvebWau7d5hZK3Ag4neJSJ3I3mROkSCyN6Np6GlggZmda2aNwE3AmmDbGmBZsLwMGEoNQ0QkS3cdjS7q8NEbzGwP8B7gfjN7IEifaWZrAdy9F7gVeADYCvzU3bcEH3EHsNjMtgOLg3URkUHlho9KVJFmFrv7amB1ifR9wJK89bXA2hL5DgJXRSmDiNQnTSirHs0sFpFYSuiBNFWjQCAisWRqGqoaBQIRiaXsRDLVCCJTIBCRWMo1DdW4IKOAAoGIxFJuQpkiQVQKBCISS6abzlWNAoGIxJKahqpHgUBEYinXV6xIEJUCgYjEkmoC1aNAICKxpAll1aNAICKxpM7i6lEgEJFYSugu1FWjQCAisaSmoepRIBCRWMpOKFMciEyBQERiScNGq0eBQERiKewjUNNQdAoEIhJLCY0aqpqoj6q80cy2mFnazNrK5JljZo+Y2dYg72fytt1uZnvNbEPwWlLqM0REilnRu1Qu0qMqgc3AR4H/N0CeXuDv3f2PZtYCPGNmD7n788H2b7j7VyOWQ0TqTG4egUJBVFGfWbwVBj4R7t4BdATLx8xsKzALeL7sTiIig9BN56rnTe0jMLN5wMXAk3nJt5rZJjO7x8wmvZnlEZH40hPKqmfQQGBmD5vZ5hKvpcP5IjMbD/wc+Ky7Hw2SvwOcBywiU2v42gD7LzezdjNr7+zsHM5Xi8golEjomcXVMmjTkLtfHfVLzKyBTBD4sbv/Iu+z9+fluRv49QDlWAmsBGhra/OoZRKReDMNH62aM940ZJkOhO8BW93960XbWvNWbyDT+SwiMigNH62eqMNHbzCzPcB7gPvN7IEgfaaZrQ2yXQ7cDHyoxDDRO83sOTPbBFwJfC5KeUSkfuimc9UTddTQamB1ifR9wJJgeT1lzpW73xzl+0Wkfummc9WjmcUiEkvZYeuKA5EpEIhILKlpqHoUCEQklgw1DVWLAoGIxFK2RqA4EJkCgYjEUnZCmQJBZAoEIhJLmlBWPQoEIhJLCgDVo0AgIrGU6yNQQIhKgUBEYkm3oa4eBQIRiaXsg2k0kyAyBQIRiSUNH60eBQIRiSU1DVWPAoGIxFKuJqBIEJUCgYjEkul5BFWjQCAisaSbzlWPAoGIxJKeR1A9CgQiEksaNVQ9CgQiEkvqI6ieqM8svtHMtphZ2szaBsi3K3g28QYza89Ln2xmD5nZ9uB9UpTyiEj9yD28XpEgqqg1gs3AR4FHh5D3Sndf5O75AWMFsM7dFwDrgnURkUGps7h6IgUCd9/q7tsifMRSYFWwvAq4Pkp5RKR+qEZQPW9WH4EDD5rZM2a2PC99hrt3AATv08t9gJktN7N2M2vv7Ow8w8UVkZEu9zyC2pZjNEgNlsHMHgbOLrHpNne/b4jfc7m77zOz6cBDZvaCuw+lOSnL3VcCKwHa2tp8OPuKyOijm85Vz6CBwN2vjvol7r4veD9gZquBS8n0K+w3s1Z37zCzVuBA1O8Skfqg4aPVc8abhsxsnJm1hMvANWQ6mQHWAMuC5WXAUGsYIlLnEho+WjVRh4/eYGZ7gPcA95vZA0H6TDNbG2SbAaw3s43AU8D97v7bYNsdwGIz2w4sDtZFRAaVGzWkSBDVoE1DA3H31cDqEun7gCXB8k7gojL7HwSuilIGEalPmlBWPZpZLCKxpHsNVY8CgYjEkjqLq0eBQERiKdtZXONyjAYKBCISS5atESgURKVAICKx9CfzJvPpK+bzjlln1boosRdp1JCISK2Ma0rxxSVvr3UxRgXVCERE6pwCgYhInVMgEBGpcwoEIiJ1ToFARKTOKRCIiNQ5BQIRkTqnQCAiUufMPX5PfTSzTuCVCnefCrxexeLUko5lZNKxjEw6FjjH3acVJ8YyEERhZu3u3lbrclSDjmVk0rGMTDqW8tQ0JCJS5xQIRETqXD0GgpW1LkAV6VhGJh3LyKRjKaPu+ghERKRQPdYIREQkjwKBiEidq6tAYGbXmtk2M9thZitqXZ7hMrNdZvacmW0ws/YgbbKZPWRm24P3SbUuZylmdo+ZHTCzzXlpZctuZl8MztM2M/vT2pS6vzLHcbuZ7Q3OywYzW5K3bUQeB4CZzTGzR8xsq5ltMbPPBOlxPC/ljiV258bMxpjZU2a2MTiWLwfpZ+68uHtdvIAk8BIwH2gENgILa12uYR7DLmBqUdqdwIpgeQXwf2pdzjJlvwK4BNg8WNmBhcH5aQLODc5bstbHMMBx3A58vkTeEXscQflagUuC5RbgxaDMcTwv5Y4lducGMGB8sNwAPAm8+0yel3qqEVwK7HD3ne7eDdwLLK1xmaphKbAqWF4FXF+7opTn7o8Ch4qSy5V9KXCvu59295eBHWTOX82VOY5yRuxxALh7h7v/MVg+BmwFZhHP81LuWMoZycfi7n48WG0IXs4ZPC/1FAhmAbvz1vcw8D+UkciBB83sGTNbHqTNcPcOyPxnAKbXrHTDV67scTxXt5rZpqDpKKyyx+Y4zGwecDGZvz5jfV6KjgVieG7MLGlmG4ADwEPufkbPSz0FAiuRFrexs5e7+yXAdcAtZnZFrQt0hsTtXH0HOA9YBHQAXwvSY3EcZjYe+DnwWXc/OlDWEmkj6nhKHEssz42797n7ImA2cKmZvWOA7JGPpZ4CwR5gTt76bGBfjcpSEXffF7wfAFaTqf7tN7NWgOD9QO1KOGzlyh6rc+Xu+4P/uGngbnLV8hF/HGbWQObC+WN3/0WQHMvzUupY4nxuANz9CPA74FrO4Hmpp0DwNLDAzM41s0bgJmBNjcs0ZGY2zsxawmXgGmAzmWNYFmRbBtxXmxJWpFzZ1wA3mVmTmZ0LLACeqkH5hiT8zxm4gcx5gRF+HGZmwPeAre7+9bxNsTsv5Y4ljufGzKaZ2cRgeSxwNfACZ/K81LqH/E3ujV9CZjTBS8BttS7PMMs+n8zIgI3AlrD8wBRgHbA9eJ9c67KWKf+/kama95D5C+YTA5UduC04T9uA62pd/kGO41+B54BNwX/K1pF+HEHZ3kemCWETsCF4LYnpeSl3LLE7N8A7gWeDMm8GvhSkn7HzoltMiIjUuXpqGhIRkRIUCERE6pwCgYhInVMgEBGpcwoEIiJ1ToFARKTOKRCIiNS5/w+YhBefKgkCkgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "tfp = (ss['C']+irf['C'])/(ss['L']+irf['L'])\n",
    "tfp = tfp/tfp[-1] - 1\n",
    "plt.plot(tfp)\n",
    "lov = (ss['N']+irf['N'])**(1/(calibration['sigma']-1))\n",
    "tfp2 = lov/lov[-1]-1 \n",
    "plt.plot(tfp2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "25b8edce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fe5efe99190>]"
      ]
     },
     "execution_count": 114,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAOHUlEQVR4nO3c34tc533H8fenUkQJSbFdybYsyV011UXVUogYhCG9CPUPJMVYvuiFDYmFcyEMNTi0wVXqf8CJoTGmxkakBpm4mEASIoyCYru5VeqVY8uoiuONSKqNFHuTCyfgCyHy7cUetevNSDu7Z1a76+f9gmHmnPOcmedhwG/NmVmnqpAkteuPVnoCkqSVZQgkqXGGQJIaZwgkqXGGQJIat36lJ7AUGzdurImJiZWehiStKSdPnvx1VW2av39NhmBiYoLJycmVnoYkrSlJfjFsv5eGJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxYwlBkj1J3k4yleTQkONJ8lR3/FSSXfOOr0vy4yQvjWM+kqTR9Q5BknXA08BeYCdwf5Kd84btBXZ0t4PAM/OOPwKc6TsXSdLijeMTwW5gqqrOVtVF4EVg/7wx+4Hna9YJ4LokmwGSbAU+B3xjDHORJC3SOEKwBTg3Z3u62zfqmCeBR4HfX+1FkhxMMplkcmZmpteEJUn/bxwhyJB9NcqYJHcD71XVyYVepKoOV9WgqgabNm1ayjwlSUOMIwTTwLY521uB8yOO+QxwT5KfM3tJ6e+SfHMMc5IkjWgcIXgN2JFke5INwH3A0XljjgIPdL8eug14v6ouVNVXqmprVU105/1nVX1+DHOSJI1ofd8nqKpLSR4GjgPrgOeq6nSSh7rjzwLHgH3AFPAB8GDf15UkjUeq5l/OX/0Gg0FNTk6u9DQkaU1JcrKqBvP3+5fFktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjTMEktQ4QyBJjRtLCJLsSfJ2kqkkh4YcT5KnuuOnkuzq9m9L8sMkZ5KcTvLIOOYjSRpd7xAkWQc8DewFdgL3J9k5b9heYEd3Owg80+2/BPxTVf0lcBvwD0POlSQto3F8ItgNTFXV2aq6CLwI7J83Zj/wfM06AVyXZHNVXaiq1wGq6nfAGWDLGOYkSRrROEKwBTg3Z3uaP/yP+YJjkkwAnwZ+NIY5SZJGNI4QZMi+WsyYJJ8Avg18qap+O/RFkoNJJpNMzszMLHmykqQPG0cIpoFtc7a3AudHHZPkY8xG4IWq+s6VXqSqDlfVoKoGmzZtGsO0JUkwnhC8BuxIsj3JBuA+4Oi8MUeBB7pfD90GvF9VF5IE+HfgTFX96xjmIklapPV9n6CqLiV5GDgOrAOeq6rTSR7qjj8LHAP2AVPAB8CD3emfAb4AvJXkjW7fv1TVsb7zkiSNJlXzL+evfoPBoCYnJ1d6GpK0piQ5WVWD+fv9y2JJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJapwhkKTGGQJJatxYQpBkT5K3k0wlOTTkeJI81R0/lWTXqOdKkpZX7xAkWQc8DewFdgL3J9k5b9heYEd3Owg8s4hzJUnLaByfCHYDU1V1tqouAi8C++eN2Q88X7NOANcl2TziuZKkZTSOEGwBzs3Znu72jTJmlHMBSHIwyWSSyZmZmd6TliTNGkcIMmRfjThmlHNnd1YdrqpBVQ02bdq0yClKkq5k/RieYxrYNmd7K3B+xDEbRjhXkrSMxvGJ4DVgR5LtSTYA9wFH5405CjzQ/XroNuD9qrow4rmSpGXU+xNBVV1K8jBwHFgHPFdVp5M81B1/FjgG7AOmgA+AB692bt85SZJGl6qhl+RXtcFgUJOTkys9DUlaU5KcrKrB/P3+ZbEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjDIEkNc4QSFLjeoUgyQ1JXk7yTnd//RXG7UnydpKpJIfm7H8iyU+SnEry3STX9ZmPJGnx+n4iOAS8WlU7gFe77Q9Jsg54GtgL7ATuT7KzO/wy8NdV9TfAT4Gv9JyPJGmR+oZgP3Cke3wEuHfImN3AVFWdraqLwIvdeVTVD6rqUjfuBLC153wkSYvUNwQ3VdUFgO7+xiFjtgDn5mxPd/vm+yLw/Z7zkSQt0vqFBiR5Bbh5yKHHRnyNDNlX817jMeAS8MJV5nEQOAhw6623jvjSkqSFLBiCqrrjSseSvJtkc1VdSLIZeG/IsGlg25ztrcD5Oc9xALgbuL2qiiuoqsPAYYDBYHDFcZKkxel7aegocKB7fAD43pAxrwE7kmxPsgG4rzuPJHuAfwbuqaoPes5FkrQEfUPwOHBnkneAO7ttktyS5BhA92Xww8Bx4Azwrao63Z3/b8AngZeTvJHk2Z7zkSQt0oKXhq6mqn4D3D5k/3lg35ztY8CxIeP+os/rS5L68y+LJalxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxhkCSGmcIJKlxvUKQ5IYkLyd5p7u//grj9iR5O8lUkkNDjn85SSXZ2Gc+kqTF6/uJ4BDwalXtAF7ttj8kyTrgaWAvsBO4P8nOOce3AXcC/9NzLpKkJegbgv3Ake7xEeDeIWN2A1NVdbaqLgIvdudd9nXgUaB6zkWStAR9Q3BTVV0A6O5vHDJmC3BuzvZ0t48k9wC/rKo3F3qhJAeTTCaZnJmZ6TltSdJl6xcakOQV4OYhhx4b8TUyZF8l+Xj3HHeN8iRVdRg4DDAYDPz0IEljsmAIquqOKx1L8m6SzVV1Iclm4L0hw6aBbXO2twLngU8B24E3k1ze/3qS3VX1q0WsQZLUQ99LQ0eBA93jA8D3hox5DdiRZHuSDcB9wNGqequqbqyqiaqaYDYYu4yAJF1bfUPwOHBnkneY/eXP4wBJbklyDKCqLgEPA8eBM8C3qup0z9eVJI3JgpeGrqaqfgPcPmT/eWDfnO1jwLEFnmuiz1wkSUvjXxZLUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1zhBIUuMMgSQ1LlW10nNYtCQzwC9Weh5LsBH49UpP4hpqbb3gmluxVtf8Z1W1af7ONRmCtSrJZFUNVnoe10pr6wXX3IqP2pq9NCRJjTMEktQ4Q3BtHV7pCVxjra0XXHMrPlJr9jsCSWqcnwgkqXGGQJIaZwjGKMkNSV5O8k53f/0Vxu1J8naSqSSHhhz/cpJKsnH5Z91P3zUneSLJT5KcSvLdJNdds8kv0gjvW5I81R0/lWTXqOeuVktdc5JtSX6Y5EyS00keufazX5o+73N3fF2SHyd56drNuqeq8jamG/A14FD3+BDw1SFj1gE/A/4c2AC8Ceycc3wbcJzZP5jbuNJrWu41A3cB67vHXx12/mq4LfS+dWP2Ad8HAtwG/GjUc1fjreeaNwO7usefBH76UV/znOP/CPwH8NJKr2fUm58Ixms/cKR7fAS4d8iY3cBUVZ2tqovAi915l30deBRYK9/i91pzVf2gqi51404AW5d3uku20PtGt/18zToBXJdk84jnrkZLXnNVXaiq1wGq6nfAGWDLtZz8EvV5n0myFfgc8I1rOem+DMF43VRVFwC6+xuHjNkCnJuzPd3tI8k9wC+r6s3lnugY9VrzPF9k9l9aq9Eoa7jSmFHXv9r0WfP/STIBfBr40finOHZ91/wks/+Q+/0yzW9ZrF/pCaw1SV4Bbh5y6LFRn2LIvkry8e457lrq3JbLcq153ms8BlwCXljc7K6ZBddwlTGjnLsa9Vnz7MHkE8C3gS9V1W/HOLflsuQ1J7kbeK+qTib57LgntpwMwSJV1R1XOpbk3csfi7uPiu8NGTbN7PcAl20FzgOfArYDbya5vP/1JLur6ldjW8ASLOOaLz/HAeBu4PbqLrKuQlddwwJjNoxw7mrUZ80k+RizEXihqr6zjPMcpz5r/nvgniT7gD8G/iTJN6vq88s43/FY6S8pPko34Ak+/MXp14aMWQ+cZfY/+pe/jPqrIeN+ztr4srjXmoE9wH8Dm1Z6LQusc8H3jdlrw3O/RPyvxbznq+3Wc80BngeeXOl1XKs1zxvzWdbQl8UrPoGP0g34U+BV4J3u/oZu/y3AsTnj9jH7K4qfAY9d4bnWSgh6rRmYYvZ66xvd7dmVXtNV1voHawAeAh7qHgd4ujv+FjBYzHu+Gm9LXTPwt8xeUjk1573dt9LrWe73ec5zrKkQ+L+YkKTG+ashSWqcIZCkxhkCSWqcIZCkxhkCSWqcIZCkxhkCSWrc/wLouA/ZRwywxQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(lov/lov[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "776d61a6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Solving CES Model for ['w', 'Ne', 'C'] to hit ['budget', 'free_entry', 'goods_mkt']\n",
      "Solving firms_inner for ['v', 'N'] to hit ['valfunc', 'accumulation']\n",
      "On iteration 0\n",
      "   max error for valfunc is 0.00E+00\n",
      "   max error for accumulation is 0.00E+00\n",
      "On iteration 0\n",
      "   max error for budget is 0.00E+00\n",
      "   max error for free_entry is 3.50E-02\n",
      "   max error for goods_mkt is 0.00E+00\n",
      "Solving firms_inner for ['v', 'N'] to hit ['valfunc', 'accumulation']\n",
      "On iteration 0\n",
      "   max error for valfunc is 6.11E+02\n",
      "   max error for accumulation is 9.88E+04\n",
      "On iteration 1\n",
      "   max error for valfunc is 3.53E+03\n",
      "   max error for accumulation is 3.64E-12\n",
      "On iteration 2\n",
      "   max error for valfunc is 6.28E-13\n",
      "   max error for accumulation is 3.64E-12\n",
      "On iteration 1\n",
      "   max error for budget is 3.89E+09\n",
      "   max error for free_entry is 1.40E+03\n",
      "   max error for goods_mkt is NAN\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/blocks/support/simple_displacement.py:432: RuntimeWarning: invalid value encountered in power\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "array must not contain infs or NaNs",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-115-a4162e725faf>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0mdfE\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m.035\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m.979\u001b[0m \u001b[0;34m**\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0mirf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m \u001b[0mirf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mces\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve_impulse_nonlinear\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mss\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m'w'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'Ne'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'C'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m'budget'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'free_entry'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'goods_mkt'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m'fE'\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mdfE\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/blocks/block.py\u001b[0m in \u001b[0;36msolve_impulse_nonlinear\u001b[0;34m(self, ss, unknowns, targets, inputs, outputs, internals, Js, options, H_U_factored, ss_initial, **kwargs)\u001b[0m\n\u001b[1;32m    204\u001b[0m                 \u001b[0;32mbreak\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    205\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 206\u001b[0;31m                 \u001b[0mU\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0mH_U_factored\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mapply\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mresults\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    207\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    208\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mf'No convergence after {opts[\"maxit\"]} backward iterations!'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/classes/jacobian_dict.py\u001b[0m in \u001b[0;36mapply\u001b[0;34m(self, x)\u001b[0m\n\u001b[1;32m    290\u001b[0m         \u001b[0;34m\"\"\"Returns -H_U^{-1} @ x\"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0mxsub\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mImpulseDict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtargets\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 292\u001b[0;31m         \u001b[0mout\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mfactored_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mH_U_factored\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mxsub\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    293\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mImpulseDict\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munpack\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0munknowns\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    294\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/mcr/res-m1wlg01/sequence-jacobian/src/sequence_jacobian/utilities/misc.py\u001b[0m in \u001b[0;36mfactored_solve\u001b[0;34m(Z, y)\u001b[0m\n\u001b[1;32m     47\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     48\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mfactored_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mZ\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 49\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mscipy\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlu_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mZ\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     50\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     51\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/apps/Anaconda/2021.05/lib/python3.8/site-packages/scipy/linalg/decomp_lu.py\u001b[0m in \u001b[0;36mlu_solve\u001b[0;34m(lu_and_piv, b, trans, overwrite_b, check_finite)\u001b[0m\n\u001b[1;32m    133\u001b[0m     \u001b[0;34m(\u001b[0m\u001b[0mlu\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpiv\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlu_and_piv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    134\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcheck_finite\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 135\u001b[0;31m         \u001b[0mb1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray_chkfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    136\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    137\u001b[0m         \u001b[0mb1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/apps/Anaconda/2021.05/lib/python3.8/site-packages/numpy/lib/function_base.py\u001b[0m in \u001b[0;36masarray_chkfinite\u001b[0;34m(a, dtype, order)\u001b[0m\n\u001b[1;32m    486\u001b[0m     \u001b[0ma\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0morder\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    487\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mchar\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtypecodes\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'AllFloat'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfinite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 488\u001b[0;31m         raise ValueError(\n\u001b[0m\u001b[1;32m    489\u001b[0m             \"array must not contain infs or NaNs\")\n\u001b[1;32m    490\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: array must not contain infs or NaNs"
     ]
    }
   ],
   "source": [
    "T = 300\n",
    "dfE = .035 * .979 ** np.arange(T)\n",
    "irf = {}\n",
    "irf = ces.solve_impulse_nonlinear(ss, ['w', 'Ne', 'C'], ['budget', 'free_entry', 'goods_mkt'], {'fE': dfE})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84bdf333",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f674f60ec70>]"
      ]
     },
     "execution_count": 104,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAd/klEQVR4nO3deXzU9b3v8dcnkz2EBEjCEghhl1XAgFaroOJeS217Kta2ahesrX303Gt7unBPa29PT0899fZ6ent7pK2ttVqX1q1qK+LFCipCUEC2EHbCkoWQAFknM9/7RwYOYhZgJvzmN/N+Ph55zMxvht+8f/nBm99857eYcw4REfGvFK8DiIhIdFTkIiI+pyIXEfE5FbmIiM+pyEVEfC7VizctKChwpaWlXry1iIhvrVmzps45V3jqdE+KvLS0lPLyci/eWkTEt8xsd1fTNbQiIuJzKnIREZ9TkYuI+JyKXETE51TkIiI+F5MiN7NrzazCzLaZ2bdjMU8RETk9URe5mQWAXwDXAZOAW8xsUrTzFRGR0xOL/chnA9ucczsAzOxxYD6wKQbzFhEBIBR2BENhgqEwHaHI/bCjIzItGJnWEXaEw45Q2BFyjnCYyO3J0zpvQ2FH2DlCYbqYdsrzp3HK71NPC97VH/no9GGMHJQTq18LEJsiLwb2nvS4Crjw1BeZ2UJgIUBJSUkM3lZEzpVw2NESDNHU3kFzW4jm9hAtwQ7agmHaOsK0dYRoDXbetnWEaQuGaQ2GTjzX1nHS4+D7X3+8gI+XcLCjs6CPF3Z7KExHKEw4QS6dMHV4XlwWuXUx7QO/cufcYmAxQFlZWYKsEpH4FQo7jrYGOdraQWNLkCOR+0daIretQZraOmhqD9HSHqKprYPm9hDN7Z23Te0dkekhWoKhs8qQYpCZFiAjNYWM1AAZaSlkRm4zUlPISg+QF0ghNZBCeiCF1ICRFkghLXKbmpJCWqqRlpLS+ThgXb8ukEJ6wEhNSSEQMAJmBFKMlMhtIIUT9/9r2kn3zUhJ4aT7J91G7puBndJ21kX9ffA1p/5OuqrM6MSiyKuAESc9Hg7sj8F8RYTOreEjrUHqm9o53NzOoWOdt/VNwchtO4eb2mk8qaCPtnZwrK2j13lnpwciP6kn7udkpFLQL4OcjNQPPp+RSk5kWmZa4ERJnyjrUx6nphjWB8Ul7xeLIl8NjDOzUcA+YAHw6RjMVyShHWvroOZIKzVH2zp/jt+P3NYebTtR3t0NK2SkpjAoJ5387HTys9MoLcimf2Ya/bPSyM1M7eJ+5DYzjX6ZqQRSVLKJIOoid851mNndwMtAAHjIObcx6mQiPtbWEeJAQyv7G1qoamhh3+EW9je0sK+hhQONrdQcaaWp/YPDFempKRTlZlCUm8HowhxmjRrIwOx0BuSkMyin87bzcRoDc9LJTvfkvHcSZ2Lyt8A59xLwUizmJeIXDc3t7KxrYtehJnbWNbP7UBN76pvZd7iF2mNt79tjwQyKcjMozs9i0rD+XD6hiKL+GQzun0FRbmakvDPpn5WqoQg5Y/rvXKQHobBjT30zFQePsrX6KDtqj7HzUDO76ppobAmeeJ0ZDMvLYuSgbOaML6R4QBbF+VkUD8hieH42Q/IySU/VgdTSN1TkInTu/7u/sZWKg0eoOHiMrdWdxb2t5hhtHeETryvOz6K0IJsbpg1l1KAcSgtyGFWQzYiB2WSkBjxcAklmKnJJOs459jW0sGFfI+/ta2R9VSMb9x+hvqn9xGuG9M9k/JBcLh4ziHGDc5kwOJexRf3IydA/GYk/+lspCa+5vYO1exoo332YNbsPs76qgcPNncMigRRj/OBc5k0sYmpxHhOH9mfc4FzystI8Ti1y+lTkknCqj7RSvuswq3fVs2b3YTYdOEIo7DCD8UW5XDVpMFOL85gSKe7MNA2JiL+pyMX3jrQGWbn9EG9sq2PFtjq21zYBkJmWwvQR+dw1ZwwXlA5gZskAbWlLQlKRi++Ewo61extYtqWGFdvqWF/VQNhBVlqA2aMGcvOsEcweNYjJw/qTFtCeIpL4VOTiC8faOlhRWcvSzTUs21LDoaZ2AinG+cPzuPvysVwytoAZJQO0i58kJRW5xK2G5nb+tuEgL204yMrth2gPhemfmcrcCUVcObGIueOLyMvWUImIilziypHWIK9srOYv6/ezorKOjrBj5KBsPvehkVw5cTBlpQM0XCJyChW5eK69I8yrm6t55t19vFZRS3soTHF+Fl/48ChuPH8Yk4f112HrIj1QkYtnthw8wpOrq3h27T7qm9opys3gMxeN5CPnD2XGiHyVt8hpUpHLOdXSHuLZtfv446o9rK9qJC1gXD1pCJ8sG85l4wp1WlWRs6Ail3Nib30zj6zczROr99LYEuS8Ibl8/8ZJzJ9ezMCcdK/jifiailz6jHOOlTvq+c2Knby6pZoUM66dPITbLi5lVukADZ2IxEhURW5m/wDcC0wEZjvnymMRSvwtHHa8srmaX762nbV7GxiUk85X547l1otKGJqX5XU8kYQT7Rb5BuDjwIMxyCI+FwyFeX7tfv7z79uprDnGiIFZ/MvHpvDJC4brfCYifSiqInfObQb0ETnJhcKO59ft42evVLKnvpnzhuTywILp3DB1KKna51ukz52zMXIzWwgsBCgpKTlXbyt9yDnHkk3V3L+kgq3Vx5g8rD+//lwZV04s0n/uIudQr0VuZkuBIV08tcg599zpvpFzbjGwGKCsrKyba4KLX7y5vY6f/K2CdXsbGF2Ywy8+PZPrpgwhRbsPipxzvRa5c27euQgi/rC3vpl/eXETL2+spjg/i/s+OY2PzyjWEIqIh7T7oZyWprYOfvnadhYv30FqivHNaybwhQ+P0peYInEg2t0PbwJ+DhQCL5rZWufcNTFJJnHBOcfz6/bz45e2cPBIKzfNKOZb157HkLxMr6OJSES0e608AzwToywSZ/bWN7Po2Q28vrWWqcV5/OLWGVwwcqDXsUTkFBpakQ8IhR2/fWMn9y/ZSorBvTdO4rMfKtV5UETilIpc3mdbzTHueXIt66oaueK8In74sSkU5+toTJF4piIXoPOw+kdW7uZfX9pMVnqA/7hlBjdOG6r9wUV8QEUuHGxs5Zt/WsfyyjrmTijkvk9Mo6i/vswU8QsVeZJbuqmae55aR3tHmB/dNIVPzy7RVriIz6jIk1QwFOanL1fw4Os7mFLcn5/fMpNRBTlexxKRs6AiT0IHGlu4+7F3WbP7MJ+9aCSLbpioA3tEfExFnmTe3F7H3Y+9S1swxM9vmcGN5w/zOpKIRElFniScc/xh5W7u/csmRhfk8OBnL2B0YT+vY4lIDKjIk0B7R5h7/7KRx97ew5XnFfG/F0wnNzPN61giEiMq8gR3uKmdO/+whlU76/nK3DHcc/UEHaEpkmBU5Alsb30zt/12FVWHW3hgwXTmTy/2OpKI9AEVeYLasK+RO363mrZgiEe/eCGzSnWyK5FEpSJPQMsra7nrD+/QPzOVx+66mHGDc72OJCJ9SEWeYF5cf4CvP/4uY4v68bs7Zuu84SJJQEWeQJ5+p4pvPLWOC0YO4De3z6K/9kwRSQpRXWjRzP7dzLaY2Xoze8bM8mOUS87QH1ft4Z6n1nHR6EE8/PnZKnGRJBLtFXNfAaY456YBW4HvRB9JztTDb+7iO0+/x5zxhTx0+yyy0/VBSySZRFXkzrklzrmOyMOVwPDoI8mZ+N0bO/n+8xu5atJgHvzsBTpnikgSinaL/GSfB/7a3ZNmttDMys2svLa2NoZvm7yeXL2Xe/+yiasnDeb/3jqTjFSVuEgy6vUzuJktBYZ08dQi59xzkdcsAjqAR7ubj3NuMbAYoKyszJ1VWjnh+XX7+dbT67lsfCE///QM0gKx/D9ZRPyk1yJ3zs3r6Xkzuw34CHClc04FfQ68sqma//7EWmaNHMiDn7lAW+IiSS6qb8XM7FrgW8Ac51xzbCJJT97cVsdXH32HycP685vby8hKV4mLJLtoP4//HyAXeMXM1prZf8Ygk3Sj4uBR7nxkDaUF2Tz8+dk6g6GIAFFukTvnxsYqiPTsYGMrt/92FdkZAX57x2zys9O9jiQicULfkPnA0dYgd/xuNUdagjx0+yyK87O8jiQicURHjsS5YCjMVx59h63VR3no9llMHpbndSQRiTPaIo9z339+I8sr6/jxTVOZM77Q6zgiEodU5HHs0bd389jbe/jynDF8atYIr+OISJxSkcep1bvq+f5zG5k7oZBvXjPB6zgiEsdU5HFof0MLd/1hDSMGZvPAghm6xqaI9EhfdsaZ1mCIOx9ZQ2swzOMLLyAvS/uKi0jPVORx5p+f3cB7+xr51efKGFukS7SJSO80tBJH/rymiqfWVPG1K8Zy1aTBXscREZ9QkceJbTXH+OfnNjB71EC+fuU4r+OIiI+oyONAazDE3Y+9Q2ZagP9YMINUnZJWRM6AxsjjwP98YRNbDh7ld3fM0lXvReSMadPPYy+s33/ioJ+5E4q8jiMiPqQi99CBxha++/R7zCjJ556rx3sdR0R8SkXuEecc//Sn9QRDjp99arou1SYiZy2q9jCzH5rZ+shFJZaY2bBYBUt0f1i5m+WVdSy6YSKlBTlexxERH4t2M/DfnXPTnHPTgReA70UfKfHtrGviRy9t5rLxhdx6YYnXcUTE56IqcufckZMe5gC6+HIvOkJh7nlyLemBFO77xDTMdB4VEYlO1LsfmtmPgM8BjcDlPbxuIbAQoKQkebdCH3x9B+/saeCBBdO1q6GIxESvW+RmttTMNnTxMx/AObfIOTcCeBS4u7v5OOcWO+fKnHNlhYXJeYGE7bXHeGBpJddNGcJHz9fXCSISG71ukTvn5p3mvB4DXgS+H1WiBOWc47tPv0dGWgo/mD9ZQyoiEjPR7rVy8klBPgpsiS5O4nqqvIq3d9bz3esnUpSrIRURiZ1ox8j/zcwmAGFgN/Dl6CMlntqjbfzopc3MLh3IzWW6ZJuIxFZURe6c+0SsgiSyH76wiZb2EP/68Smk6Go/IhJjOpywjy2rqOH5dfv5yuVjdKEIEekTKvI+1BoM8b3nNjCmMIe75o7xOo6IJCidxrYP/Xr5DvbWt/DoFy8kIzXgdRwRSVDaIu8jBxpb+MWy7Vw7eQiXjC3wOo6IJDAVeR/58UtbCDvHohsmeh1FRBKcirwPrNpZz/Pr9nPnZaMZMTDb6zgikuBU5DEWCjvufX4jw/IyuWvuWK/jiEgSUJHH2OOr97DpwBG+c/1EstL1BaeI9D0VeQwdbQ1y/5KtzB41kI9MG+p1HBFJEiryGPrV8p3UN7Wz6PqJOimWiJwzKvIYqTnayq+X7+CGaUM5f0S+13FEJImoyGPk569uo70jzDeunuB1FBFJMiryGNhZ18QfV+3hltkljNKFlEXkHFORx8BPl1SQnprC167U7oYicu6pyKO0vqqBF9cf4IuXjtYFI0TEEyryKP3kb1sYmJPOly4d5XUUEUlSMSlyM/uGmTkzS6qzQ63ccYg3th3iq5ePJTczzes4IpKkoi5yMxsBXAXsiT6OvzywtJLC3AxuvbDE6ygiksRisUX+M+CfABeDefnGqp31vLXjEF+eM4bMNB2KLyLeiarIzeyjwD7n3LrTeO1CMys3s/La2tpo3jYuPPDqVgr6aWtcRLzX6xWCzGwpMKSLpxYB3wWuPp03cs4tBhYDlJWV+XrrffWuet7Ydoj/ccNEbY2LiOd6LXLn3LyuppvZVGAUsC5yXpHhwDtmNts5dzCmKePMA0srKeiXzq0XjvQ6iojI2V+z0zn3HlB0/LGZ7QLKnHN1McgVt8p31bNiWx2LdJpaEYkT2o/8DD3waiWDctK59SKNjYtIfIhZkTvnShN9a3x9VQPLK+v40mWjyU4/6w8zIiIxpS3yM/Dg6zvIzUjVnioiEldU5Kdpz6Fm/vreAW69aKSO4hSRuKIiP02/XrGDQIpxxyWlXkcREXkfFflpqG9q58nyvdw0o5jB/XWGQxGJLyry0/D7t3bRGgyz8LLRXkcREfkAFXkvWtpDPPzmLuZNLGJsUa7XcUREPkBF3os/rdnL4eYgCy8b43UUEZEuqch7EAo7frV8JzNK8plVOsDrOCIiXVKR92Dp5mr21DfzpUtHEzmfjIhI3FGR9+D3b+1iWF4mV08a7HUUEZFuqci7UVl9lDe2HeLWi0aSGtCvSUTilxqqGw+/tYv01BQWzBrhdRQRkR6pyLtwpDXI0+/s48ZpwxjUL8PrOCIiPVKRd+FP5VU0t4e4/eJSr6OIiPRKRX6KcNjxyMrdzCjJZ+rwPK/jiIj0KtqLL99rZvvMbG3k5/pYBfPK65W17Kxr0ta4iPhGLK6O8DPn3E9jMJ+48PCbuyjol8F1U4Z6HUVE5LRoaOUkew4189rWWj59YQnpqfrViIg/xKKt7jaz9Wb2kJn5+jj2J8r3YMAts7XLoYj4R69FbmZLzWxDFz/zgV8CY4DpwAHg/h7ms9DMys2svLa2Nlb5Y6YjFOap8irmTihiaF6W13FERE5br2Pkzrl5pzMjM/sV8EIP81kMLAYoKytzpxvwXFlWUUvN0TYdACQivhPtXisnfyN4E7AhujjeeWL1HgpzM7j8vCKvo4iInJFo91q5z8ymAw7YBdwZbSAvHGxs5f9tqeHOOWNI03lVRMRnoipy59xnYxXES39as5ewg5vLNKwiIv6T9Juf4bDjifK9fGj0IEoLcryOIyJyxpK+yN/cfoi99S0s0C6HIuJTSV/kj6/eQ15WGtdMHuJ1FBGRs5LURV7f1M6SjdXcNKOYzLSA13FERM5KUhf5X9btpz0U5lP6klNEfCypi/zpd/dx3pBcJg3r73UUEZGzlrRFvr32GOv2NvDxmcVeRxERiUrSFvmz7+4jxWD+dBW5iPhbUhZ5OOx45t19XDK2gMH9M72OIyISlaQs8vLdh6k63MJNM7Q1LiL+l5RF/sy7VWSlBbTvuIgkhKQr8tZgiBfWH+CayYPJyYjFle5ERLyVdEW+bEsNR1s7uGnmcK+jiIjERNIV+dPv7qMwN4NLxgzyOoqISEwkVZEfbmrntYoa5p8/jFSdd1xEEkRStdlfNxwkGHJ8THuriEgCibrIzexrZlZhZhvN7L5YhOorL763n9EFOUzWIfkikkCi2m3DzC4H5gPTnHNtZha3F7ysO9bGW9sP8dXLx2JmXscREYmZaLfI7wL+zTnXBuCcq4k+Ut/424aDhB3cMG1o7y8WEfGRaIt8PHCpmb1tZn83s1ndvdDMFppZuZmV19bWRvm2Z+7F9QcYU5jDhMG55/y9RUT6Uq9DK2a2FOjqEMhFkT8/ALgImAU8aWajnXPu1Bc75xYDiwHKyso+8Hxfqjnayts7D3H3FeM0rCIiCafXInfOzevuOTO7C3g6UtyrzCwMFADnfpO7By9HhlU+omEVEUlA0Q6tPAtcAWBm44F0oC7KecbcC+sPMK6oH+M1rCIiCSjaIn8IGG1mG4DHgdu6GlbxUs2RVlbtqteXnCKSsKLa/dA51w58JkZZ+sRfNxzEObhhqopcRBJTwh/Z+eL6A4wf3I9xGlYRkQSV0EV+sLGV1bvruWHqMK+jiIj0mYQu8pc3RoZVpukCEiKSuBK6yJdsOsiYwhzGFmlYRUQSV8IWeWNzkLd31HPVJG2Ni0hiS9giX1ZRQ0fYcfXkwV5HERHpUwlb5Es2HaQwN4Ppw/O9jiIi0qcSsshbgyH+XlHLvImDSUnRuVVEJLElZJG/tf0QTe0hDauISFJIyCJfsqmanPQAF+sCyyKSBBKuyMNhx9LN1cydUERGasDrOCIifS7hinxtVQO1R9u4apKGVUQkOSRckS/ZWE1qinH5hLi9fKiISEwlXJG/sukgF40eRF52mtdRRETOiYQq8u21x9he26RhFRFJKlGdj9zMngAmRB7mAw3OuelRZjprSzdVAzBPRS4iSSTaC0vcfPy+md0PNEadKArLKmo4b0guxflZXsYQETmnYjK0Yp2Xpv8U8MdYzO9sHGkNUr7rMJefpy85RSS5xGqM/FKg2jlXGaP5nbE3KuvoCDvtrSIiSafXoRUzWwp0dS7YRc655yL3b6GXrXEzWwgsBCgpKTnDmL17raKW3MxUZpbkx3zeIiLxrNcid87N6+l5M0sFPg5c0Mt8FgOLAcrKytwZZOyVc45lFTVcNq6Q1EBC7YgjItKrWLTePGCLc64qBvM6K5sOHKHmaBtzJxR6FUFExDOxKPIFePglJ3QOqwDMUZGLSBKKavdDAOfc7THIEZVlW2qYWpxHUW6m11FERM453w8oNzYHeWfPYS7X1riIJCnfF/nrlbWEHczRbocikqR8X+TLKmrIz05j+oh8r6OIiHjC10UeDjv+XlHLnPGFBHRtThFJUr4u8vf2NXKoqV1Hc4pIUvN1kb9WUYsZXDZeX3SKSPLydZGv2FbL1OI8Buakex1FRMQzvi3yY20dvLungUvGFngdRUTEU74t8lU7D9ERdlyqIheRJOfbIl9eWUdGagozRw7wOoqIiKd8W+RvbKtj9qiBZKYFvI4iIuIpXxZ5zZFWtlYf48MaVhER8WeRr9hWB6AvOkVE8HGRD8xJZ9LQ/l5HERHxnO+K3DnHiso6Lh4ziBQdli8i4r8i31ZzjJqjbRofFxGJiKrIzWy6ma00s7VmVm5ms2MVrDvHx8c/PE5FLiIC0W+R3wf8wDk3Hfhe5HGfWlFZR+mgbIYPyO7rtxIR8YVoi9wBx79xzAP2Rzm/HgVDYVbuOKS9VUREThLtNTv/EXjZzH5K538KF3f3QjNbCCwEKCkpOas3W7e3gab2EJdqWEVE5IRei9zMlgJDunhqEXAl8N+cc382s08BvwHmdTUf59xiYDFAWVmZO5uwyyvrMIMPjVaRi4gc12uRO+e6LGYAM/s98PXIw6eAX8coV5eG5WfyDxcMJy87rS/fRkTEV6IdWtkPzAFeA64AKqMN1JObZ5Vw86yzG5YREUlU0Rb5l4AHzCwVaCUyBi4iIudOVEXunFsBXBCjLCIichZ8d2SniIi8n4pcRMTnVOQiIj6nIhcR8TkVuYiIz6nIRUR8zpw7q6Plo3tTs1pg91n+8QKgLoZxvKRliU9alvikZYGRzrnCUyd6UuTRMLNy51yZ1zliQcsSn7Qs8UnL0j0NrYiI+JyKXETE5/xY5Iu9DhBDWpb4pGWJT1qWbvhujFxERN7Pj1vkIiJyEhW5iIjP+arIzexaM6sws21m9m2v85wpM9tlZu+Z2VozK49MG2hmr5hZZeR2gNc5u2JmD5lZjZltOGlat9nN7DuR9VRhZtd4k/qDulmOe81sX2S9rDWz6096Li6XA8DMRpjZMjPbbGYbzezrkel+XC/dLYvv1o2ZZZrZKjNbF1mWH0Sm9916cc754gcIANuB0UA6sA6Y5HWuM1yGXUDBKdPuA74duf9t4Cde5+wm+2XATGBDb9mBSZH1kwGMiqy3gNfL0MNy3At8o4vXxu1yRPINBWZG7ucCWyOZ/bheulsW360bwIB+kftpwNvARX25Xvy0RT4b2Oac2+GcawceB+Z7nCkW5gMPR+4/DHzMuyjdc869DtSfMrm77POBx51zbc65ncA2Otef57pZju7E7XIAOOcOOOfeidw/CmwGivHneuluWboTz8vinHPHIg/TIj+OPlwvfiryYmDvSY+r6HlFxyMHLDGzNWZ2/LJ4g51zB6DzLzNQ5Fm6M9dddj+uq7vNbH1k6OX4R17fLIeZlQIz6Nz68/V6OWVZwIfrxswCZrYWqAFecc716XrxU5FbF9P8tu/kJc65mcB1wFfN7DKvA/URv62rXwJjgOnAAeD+yHRfLIeZ9QP+DPyjc+5ITy/tYlpcLU8Xy+LLdeOcCznnpgPDgdlmNqWHl0e9LH4q8ipgxEmPhwP7PcpyVpxz+yO3NcAzdH58qjazoQCR2xrvEp6x7rL7al0556oj//DCwK/4r4+1cb8cZpZGZ/E96px7OjLZl+ulq2Xx87oBcM41AK8B19KH68VPRb4aGGdmo8wsHVgAPO9xptNmZjlmlnv8PnA1sIHOZbgt8rLbgOe8SXhWusv+PLDAzDLMbBQwDljlQb7TcvwfV8RNdK4XiPPlMDMDfgNsds79r5Oe8t166W5Z/LhuzKzQzPIj97OAecAW+nK9eP0N7xl+G3w9nd9mbwcWeZ3nDLOPpvOb6XXAxuP5gUHAq0Bl5Hag11m7yf9HOj/aBuncgvhCT9mBRZH1VAFc53X+XpbjEeA9YH3kH9XQeF+OSLYP0/kRfD2wNvJzvU/XS3fL4rt1A0wD3o1k3gB8LzK9z9aLDtEXEfE5Pw2tiIhIF1TkIiI+pyIXEfE5FbmIiM+pyEVEfE5FLiLicypyERGf+/9lueNR0MWnDwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(100*((irf['N'] + ss['N'])/ss['N']-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42e2f2ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "tfp = (irf['C'] + ss['C'])/(irf['L'] + ss['L'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d01182b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjgAAAEYCAYAAABRMYxdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABYX0lEQVR4nO3dd3gc1dX48e9Rt9xtuTe5V4wxxsbG9GYI4AAJJUCAkJeQN71D4EdIfwMJEBIS4hDY0DHGphfbFBuDDS7IvclykWxZ3eqypNX5/TEjs5ZVVtKutuh8nmcf7ZSdPSvdvTpz5869oqoYY4wxxkSTmFAHYIwxxhgTaJbgGGOMMSbqWIJjjDHGmKhjCY4xxhhjoo4lOMYYY4yJOpbgGGOMMSbqWIITIiLybRHJEZEyEekbohhSRURFJK6J7feJyDPu8+FurLEtHPMWEVkVjHiNaUxL5dgY0zlFXYIjIvtEpFpEUhqsT3MrwdQQheYbSzzwIHCRqnZT1YIAHTdoyYWqHnBj9Qbj+MHg/r3HhDoO84VI+H4aE2zuyWL9o05EKn2Wb3BPLmsa7Pdz97UfikiVuy5fRBaLyKBQf6ZwFHUJjmsvcH39goicBHQJXTgnGAAkAVtDHYgxIRCw76e12phI5J4sdlPVbsAB4HKfdc+6u73ou5+q3u9ziO+6rx0H9AIe6thPEBmiNcF5Gvi6z/LNwFO+O4jIl0TkcxEpEZFMEbnPZ1uSiDwjIgUickRE1orIAHfbLSKSISKlIrJXRG5oLAARSRSRh0XkkPt42F03Dtjp7nZERN5v4vWni8gn7vtvFJFzfLadEIOITAQeA2a7mf2Rlj6nj2+4MWaLyE+aiOe4ywAt/R5E5M8iUuRuu8Rn/Yci8jv3s5WJyOsi0ldEnnVjXOt7Fi8iE0RkmYgUishOEbnGZ5tHRB4VkTfdOD4VkdHutpXubhvd97lWRFJE5A33d1ooIh+JSLR+B8JZs9/PFr6b9eXwNhE5AJzw/RGRq92WoiluGfmdz7ZzRCTLZ3mfiNwlItvc8vqkiCQF+PMaExSqWgi8DEwJdSzhKFor9zVADxGZKE6fkWuBZxrsU45TyfYCvgR8W0S+7G67GegJDAP6AncAlSLSFXgEuERVuwNzgLQmYrgbOB2YBpwMzATuUdVdwGR3n16qel7DF4rIEOBN4HdAH+CnwMsi0q+pGFR1uxvnajfb7+XH56x3LjAWuAi4U0QuaOIz1cfX0u9hFk4SlwLcD/xHRMRn+3XATcAQYDSwGnjS/azbgV/5vM8y4DmgP85Z/z9EZLLPsa4Hfg30BtKB3wOo6lnu9pPd38eLwE+ALKAfTivaLwGbq6TjtfT99KfMng1MBC72XSkitwJ/Ai5Q1S1+xnODe5zROGfE97TmwxgTKuJc6r0a+DzUsYSjaE1w4IuzxAuBHcBB342q+qGqblbVOlXdBDyPU2kC1OAkNmNU1auq61W1xN1WB0wRkS6qmq2qTV1mugH4jarmqmoezj/hm/yM/UbgLVV9y41vGbAOuLSVMbT0Oev9WlXLVXUzTqJx/QkHOlFzMexX1X+7/XX+CwzCSSjqPamqe1S1GHgb2KOqy1W1FngJOMXd7zJgn6o+qaq1qroB52zlKz7HWqyqn7mvfRYnoWxKjRvLCFWtUdWP1CZjC5Umv59+ltn73DJb6bPuh8DPgHNUNb0VsfxdVTPds+Hf41/5NybYrnFbm+sfg322PeK20m8EsoEfhyTCMBftCc7XgFtocHkKQERmicgHIpInIsU4rR8pPq99F3jBvXRzv4jEq2o5ztnmHUC2e2lkQhPvPxjY77O8313njxHAV30LNzAXGNTKGFr6nPUyWxOnHzEc9tm3wn3azWd7js/zykaW6/cdAcxq8Hu4ARjY2HsBFQ3ep6EHcFp5lrqX1+5sZl8TXE1+P9tQZuv9DHhUVbMa2dacVpV/YzrIQlXt5fM45LPt++66Iap6g3sSbRqI2gRHVffjdGa8FFjcyC7PAa8Bw1S1J07/FXFfW6Oqv1bVSTiXXy7D7TOgqu+q6oU4LQE7gH83EcIhnH/Q9Ya76/yRCTzdoHB3VdX/ayGGxlojmvycPoa1Ns5W/B7aIxNY0eD30E1Vv92Wg6lqqar+RFVHAZcDPxaR8wMasfFLC99Pf8psY2X9IuAeEbnaZ105kOyzPJATtbr8G2PCX9QmOK7bgPPcFoeGugOFqlolIjNxziYBEJFzReQkt39ACc6lDa+IDBCRK9y+IUeBMqCp26afx6ls+7nXSe/lxH5ATXkGuFxELhaRWHE6PZ8jIkNbiCEHGCoiCf58Th//T0SS3b4ttwIvNhdcK38P7fEGME5EbhKRePdxmjgdqv2RA4yqXxCRy0RkjNsfqMSNOWJue49CTX0//SmzjdkKzAMeFZEr3HVpwKUi0kdEBuJcxmroO+53qw9Ov6xmy78xJjJEdYLj9vNY18Tm/wV+IyKlOMnHQp9tA4FFOP8EtwMrcJKOGJyOqoeAQpx+Af/bxPF/h9NvZhOwGdjgrvMn7kxgPk5lm4fTkvEz9/2bi+F9nEr+sIjk+/E5663AuXTzHvBnVV3aQoit+T20maqW4pyVX+e+12GcDqSJfh7iPuC/7uWta3A6Ui/HSchWA/9Q1Q8DHLbxUzPfT3/KbFPH3IjT4vpvce7eexqnn8I+YCmNJy/Pudsy3Idf31NjTHgT62NpjOmsRGQf8E1VXR7qWIwxgRXVLTjGGGOM6ZxaTHBE5Ay3rwUicqOIPCgiI1p6nTHRTkSeEJFcEWl0vBVxPCIi6SKySUSm+2ybJ87Ahel2N5cxxgSePy04/wQqRORk4Oc4t1GecNu1MZ2QB6dTa1Muwen3Mxa4Hee7hNt5/VF3+yTgehGZFNRITaNUNdUuTxkTnfxJcGrdwdDmA39V1b/i3OVgTKemqitxOlk3ZT7wlDrWAL3EmRRvJpCuqhmqWg284O5rjDEmQPyZqK5URO7CGV33LPfsMz64YbVPSkqKpqamhjoME+bWr1+fr6r9gvgWQzh+ELksd11j62c1dRARuR2nBYiuXbueOmFCk+M6GnNMB5TvgLJ62/jL37LtT4JzLc44FLep6mERGY4zImy7icg84K9ALPB4/UB2PtvF3X4pzii1t7jD9TcrNTWVdeuaujvcGIeI7G95r/a9RSPrtJn1jVLVBcACgBkzZqiVbeOPDijfAWX1tvGXv2Xbn0tUpTiXpj4SZybsaTiD2LWLn/0QGu3DYEyEyOL4UXKH4ozn09R6Y6KGdaQ3oeZPgrMSSBRnhuv3cEa69QTgvf3ph9BUH4ZmFRQUkJaWBoDX68Xj8bBp0yYAampq8Hg8bNni3PhSVVWFx+Nh+/btAFRUVODxeNi5cycAZWVleDwe0tOdufuKi4vxeDxkZGQAUFRUhMfjYd++fQDk5+fj8XjIzHSuQOTm5uLxeDh40JlL8PDhw3g8Hg4fdqZQOnjwIB6Ph9zcXAAyMzPxeDzk5zvj9O3btw+Px0NRUREAGRkZeDweiouLAUhPT8fj8VBWVgbAzp078Xg8VFQ4U0Bt374dj8dDVVUVAFu2bMHj8VBTUwPApk2b8Hg8eL3OgL5paWl4PJ5jv8v169fz1FNf9Clfu3Ytzz777LHlNWvW8PzzX+S7n3zyCQsXfjEu26pVq1i0aNGx5RUrVrB48Rcj83/wwQe8+uqrx5aXL1/O66+/fmx56dKlvPnmm8eW33nnHd55551jy2+++SZLly5FVfHWKa++9hpLly2jqsZLZbWXlxcvYdny947tv3jxYlasWEEHeQ34uns31elAsapmA2uBsSIy0h11+jp3X2OignWkN/5QVerqnLq7xlvH0VovVTVfPNo7Tp8/l6hEVStE5Dbgb6p6v4iktetdHf70Q2iqD0P2CUH69FMYMmRIAMIz/vDWKdW1dezJK6Oksoa9+eUUllTxatpBKqu9HN5fSGVxGX9ZupPKai9Hs3Koqyxl+bPrqa5VuhUeJKa2kmcfW021t44BZVnEeKt5cOeH1HjrGF29H6mr487P3sWryjTZhyp8/6O3qVPltNj91GoMt7//FgBz4jM5qnHc/l41AHPjD6LxhVx4QeA/u4g8D5wDpIhIFvAr3P5pqvoY8BbO5dV0nEust7rbakXkuzgTusYCTzQ3I7wxEejYCSyAiNSfwG5r6gX1J6bTpk3D6/Xy9NNPM336dKZOnUpNTQ3PPvssM2bMYMqUKVRVVfHCCy8wa9YsJk6cSEVFBQsXLmT27NmMHz+esrIyFi1axNy5cxkzZgzFxcUsWbKEs846i1GjRlFUVMSrr77KOeecQ2pqKvn5+bzxxhucf/75DBs2jNzcXN566y0uvPBChgwZwuHDh3nnnXeYN28eAwcO5ODBgyxbtoxLL72U/v37k5mZyXvvvcdll11GSkoK+/bt48MPP2T+/Pn07t2bjIwMVq5cyZVXXknPnj1JT09n1apVfOUrX6Fbt27s3LmT1atXc80115CcnMz27dv59NNPue6660hKSmLLli2sW7eOG264gfj4eDZt2sSGDRu46aabiI2NJS0tjbS0NG655RbAOTHdunUrX//61wHnxHTXrl3ccMMNqCqrPlnNnoy9zL7wCkqqatixcR1FudkMm3E+FdVecndtpKqkgNjRs6ms8VKTtZW6yhIK+51CbV0dXQt3ElNTSVaPyVR7lYGlu4mpq2ZHwjiqa+sYU70H0TrW16VSW6ecEuPU2+u8w1GF02L3UUsMa2uchuw58fs4qnGsrx0KwNz4vZRrAs/e+w2SE+JYvHgxffv25eyzz25VIfQrwRGR2TizON/mrott1bs0cdxG1jVM1/zuq9Cwn8K0adMAiI2NPfZHB4iPjz9uOSkp6bjl5OTk45a7det23HLPnj2PW+7du/dxyykpKcct9+/f/7jlgQMHHrc8ZMiQ45aHDRt23HJqaupxy6NGjWLUqGPTKzFmzBjGjBlzbHn8+PGMHz/+2PLEiROZOPGLqZumTJnClClTji1PnTqVqVOnHlueNGUqfYePY8OBIvJKj5JXk0LBoNP55ZLN5JUe5UhFDcWVQ3nkD8spqaylssYL9IKtvq0ifSE9zX0eC/RDMtPpEh9LckJvkuJTSMopIz42hoT4VOKTYoiPEXokxOPtNhmJiWFKXAzxsUJCbArxsTGMjBXiYoQYGQEC00SIEYiR0YgIZwgIQoyMIyZGON9dhgl0if+iofKqq64iUFT1+ha2K/CdJra9hZMAGRON/OpIbyemgVFR7bRYf5yeT37ZUfbtzqe4oJwfL0wjv6ya+MIMulYV8OBvl1FSVcNYDjMotoT7tq4EYErcYfrFlPHgno0AnBSXSx+p4OODGXRJiOWk2HJ6Sg0Z+U69PaymjiSF+NgYkhNi6FYTS1xdPCcN6kl8bAxdcrOJpY6RI4YSFyPUZTo3m05LHUWMQPXeIiQmljkjxyIClelFxMQnct6o8YhA6Y5CYhKTiY9t31jELU7VICJnAT8FPlbVP4nIKOCHqvr9dr2xkzTdp6oXu8t3AajqH332+Rfwoao+7y7vBM5xm/mbZB0xm1daVcPu3DL25JaRWVRJVlEFWYWVZBZVcLikisaKRN+uCfTrnkiv5Hh6dnEePZLcn/XLXeLolhhPckIsSfGxJCfE0iU+li4JsSTGxeD0GQ8fIrJeVWeEOo7WsLJt/BXK8i0iXwUuVtVvuss3ATNV9XtNvcbKdtNUlYNHKtmdW8a+/HIyC516u77+Lq2qPeE1CbEx9OueSEq3BHolJ9CjSzw9kuKOr7OT4umeFEfXxLjj6usu7vP2JhjB4m/ZbrEFxx3rY6XPcgbQruTGdawfAnAQpx9Cw1mDXwO+6zZvzuKLPgzGD6rKgcIK0jKPsPVQCTsPl7I7p5RDxVXH9hGBgT2SGNY7mdmj+zK0dzJDeiXRv3sS/bon0q97In26JoRtQTfGhCXrSN9GVTVeth4qJi2zmO3ZJezOLSM9p5Tyau+xfbrExzK0dxeG9UnmtNTeDOudzICeSfTvnkhKN6fe7pEUF3YnlR3Nn0tUQdFUPwQRucPd3mQfBtO4Gm8dm7KO8El6ARsOFLExq5jCcqcvSkJcDKP7dWPmyD6MHdCdcQO6M6Z/N4b06kJCnCUvxpiA8ucE1gAFZUf5ZE8Bn+4tIC3zCDuyS6mtc5rR+3VPZNyAbnx1xjDGDujGuAHdGZnSlb5dEzp98uKPkCU40Hg/BDexqX/eZB8G48gsrOC97TmsSs9nTUYhZUdrEYEx/bpx/oT+TBvei2nDejF+QHfirBXGGNMBrCN902q9dXy2r5APd+axanc+27JLAOiWGMfUoT25/axRnDysF6cM60X/HkkhjjayhTTBMW2zL7+ct7Zk8/bmw2w+6NwuPqJvMvOnDWbumBRmj+5Lr+SEEEdpjOnMrCP9F2q8daxKz+edzYdZtj2HwvJqEmJjOHVEb3560TjOGJPCSUN62klogLWY4LhNjN8DUn33V9UrgheWaaiiupY3N2Xz4tpM1u13xsQ5eVgv7rpkAvOmDGRE364hjtAYY4yvPXllLFybycsbDpJfdpRuiXGcN6E/l0wZyNnj+5GcYG0MweTPb/cV4D/A60BdUKMxJ8gsrOA/q/ayaH0WZUdrGZXSlTsvmcDlJw9mSK8uoQ7PGGOMj7o65cNduSxYmcGajEJiY4TzJvTnq6cO5axx/UiKD8QoK8Yf/iQ4Var6SNAjMcfZdqiEf63cwxubshHg8pMHc/3M4ZyW2ts6lxljTJip9dax5PODLFiZwe7cMgb3TOIX8yZw9alD6N/d+tKEgj8Jzl9F5FfAUuBo/Up/Jr00rZdZWMGfl+7k1bRDdE2I5RtnpHLrGSMZbK01xhgTdlSVd7Yc5oF3d5KRX86Egd156NqTuWzqYBteI8T8SXBOAm4CzuOLS1TqLpsAKa2q4eHlu3lq9T5iY4T/PWc03zprND2T40MdmjHGmEZszirm/726hbTMI4zp340FN53KhZMGWCt7mPAnwbkSGOVOiGkCTFV5d2sO9722lZzSKq6dMYwfXjCOgT2tSdMYY8JRaVUNf1m6i6dW76Nvt0Tuv3oqV00fYndBhRl/EpyNQC8gN7ihdD6F5dXc+fImlm7LYcLA7vzzxumcMrx3qMMyxhjThE/S8/nJSxs5XFLFjbNG8NOLx9Ozi7W0hyN/EpwBwA4RWcvxfXDsNvF2+Dg9nx+9mMaRihruvGQCt80daddrjTEmTNV463hw2S4eW7GHkSldefnbc5huJ6RhzZ8E51dBj6IT8dYpDy3bxaMfpjMqpStP3noakwf3DHVYxhhjmpBdXMkdz2xgY+YRrjttGPdePsnGsIkA/ky2uUJEBgCnuas+U1W7XNUGpVU1/PCFNN7bkcs1M4Zy3xWT7UtijDFh7PMDRdz+9Hoqjtby6Nem86Wpg0IdkvFTi9dEROQa4DPgq8A1wKci8pVgBxZtsooquPqfn/Dhrjx+O38y93/lZEtuooCIzBORnSKSLiJ3NrL9ZyKS5j62iIhXRPq42/aJyGZ327qOj94Y05zXNx7i2gVr6BIfy5LvnGHJTYTx5z/s3cBp9a02ItIPWA4sCmZg0SQ9t4yb/vMpZUdr+e+tM5k7NiXUIZkAEJFY4FHgQiALWCsir6nqtvp9VPUB4AF3/8uBH6lqoc9hzlXV/A4M2xjjh2c/3c89r2zhtBF9eOymU+nT1eb3izT+JDgxDS5JFeBHy49xbDlYzM1PfIYIvHj7bCYN7hHqkEzgzATSVTUDQEReAOYD25rY/3rg+Q6KzRjTRv9asYc/vr2D8yb05x83TLfpFSKUP4nKOyLyrojcIiK3AG9iM8T6ZXt2CV/79xqS4mN56Y45ltxEnyFAps9ylrvuBCKSDMwDXvZZrcBSEVkvIrc39SYicruIrBORdXl5eQEI2xjTlAUrneTm8pMH86+bTrXkJoI124IjznCMj+B0MJ4LCLBAVZd0QGwRLSPPuSzVNTGOF24/nWF9kkMdkgm8xoYr1Sb2vRz4uMHlqTNU9ZCI9AeWicgOVV15wgFVFwALAGbMmNHU8Y0x7fT8Zwf4w1s7uGzqIB6+dhqxMTYicSRrNsFRVRWRV1T1VGBxB8UU8Q4eqeTGxz9FFZ755ixLbqJXFjDMZ3kocKiJfa+jweUpVT3k/swVkSU4l7xOSHCMMcH3xqZD/HLJZs4d348Hr7HkJhr4c4lqjYic1vJuBqD8aC23edZSWlXLU7fNZHS/bqEOyQTPWmCsiIwUkQScJOa1hjuJSE/gbOBVn3VdRaR7/XPgImBLh0RtjDnO5weK+PHCjcwY0Zt/3HAqCXHWzTQa+NPJ+FzgWyKyHyjHaZZXVZ3a1jcVkQdwmuyrgT3Arap6pJH99gGlgBeoVdUZbX3PjlBXp/zghTR255bx5C02gF+0U9VaEfku8C4QCzyhqltF5A53+2PurlcCS1W13OflA4Al7qR8ccBzqvpOx0VvjAFnEL/bn17PwB5JLLhpBl0SrM9NtGgywRGRkaq6F7gkCO+7DLjL/QfxJ+Au4BdN7Bsxt9He/+5Olm/P4ddXTOascf1CHY7pAKr6Fg063fskNvXLHsDTYF0GcHKQwzPGNKOy2sv/PLWOymovz35zFr3tVvCo0lwLziLgVJyz0vMD+aaqutRncQ0Q8QMHLtuWw2Mr9vC1WcP5+uwRoQ7HGGNMC37zxja2HCzhPzfPYNyA7qEOxwRYcwlOjIj8ChgnIj9uuFFVHwxQDN8AXmxiW/1ttAr8y72bpFHubba3AwwfPjxAofnn0JFKfrZoI5MH9+BXl0/CvexgjDEmTL25KZvnPzvAHWeP5vyJA0IdjgmC5hKc64Avu/u0OrUVkeXAwEY23a2qr7r73A3UAs82cRi/bqOF0N1KW+ut44cvpFFTW8ffvzadxDi7fmuMMeEss7CCOxdvYtqwXvzkonGhDscESZMJjqruBP4kIptU9e3WHlhVL2huu4jcDFwGnK+qjSYkkXAb7T8+3MNn+wp5+NppjEzpGupwjDHGNKOuTvnRi2mg8LfrTyE+1u6YilYt/mXbkty0RETm4XQqvkJVK5rYJ+xvo915uJS/vb+by08ezJdPaXQAW2OMMWHk6TX7Wbe/iF9dMdnGKItyoUpd/45z2WuZO5PyYwAiMlhE6u9IGQCsEpGNOLOZvxlOt9F665Sfv7yJ7knx3Hf5pFCHY4wxpgUHj1Ry/zs7OHNsCldPt5PSaOfPODgBp6pjmlh/CLjUfR7Wt9E+sWovGzOP8Mj1p9C3W2KowzHGGNMMVeXuJZtR4A9XnmQ3g3QCfiU4IjIHSPXdX1WfClJMYS+rqIK/LNvJBRMHcPnUQaEOxxhjTAve3JzNhzvzuPeySXZpqpNoMcERkaeB0UAazojC4Ny+3WkTnD++vQOAX8+fbGcBxhgT5qpqvPzxrR1MHNSDm+ekhjoc00H8acGZAUxq6k6nzubTjALe3JTNDy8Yy5BeXUIdjjHGmBb8Z9VeDh6p5IGvTrVJNDsRfzoZb6Hx8Ww6HW+d8uvXtzG4ZxLfOmt0qMMxxhjTgtySKh79IJ2LJw9gzuiUUIdjOpA/LTgpwDYR+Qw4Wr9SVa8IWlRh6qV1mWzLLuFv159iE7IZY0wEeODdndR46/jlpRNDHYrpYP4kOPcFO4hIUFXj5eHluzlleC8us47FxhgT9tJzy3h5QxbfOGMkI/raQKydjT8D/a0AduCMW9Md2O6u61Se/+wAh0uq+OlF461jsTHGNEFEvioiW0WkTkRmhDKWR97bTVJ8LN8+x7oUdEYtJjgicg3OQHtfBa4BPhWRiJ/9uzUqq708+sEeZo3sw5zRfUMdjjHGhLMtwFWEeFqdXTmlvL7pELfMSbWxyjopfzoZ3w2cpqo3q+rXceaD+n/BDSu8PL1mH/llR/mJtd6YBkRknojsFJF0Ebmzke3niEixO2J3mojc6+9rjYlEqrrdncswpP66fDddE+L4nzNHhToUEyL+9MGJUdVcn+UCQjfFQ4erqK7lsRUZnDk2hZkj+4Q6HBNGRCQWeBS4EMgC1orIa6q6rcGuH6nqZW18rTFRS0RuB24HGD58eMCOm55bxltbsvnOOWPo3TUhYMc1kcWfBOcdEXkXeN5dvhZ4q5n9o8rCtZkUllfzwwvGhjoUE35mAunutCKIyAvAfMCfJKU9rzUmpERkOY0PH3K3qr7q73FUdQGwAGDGjBkBG2vt8Y8ySIiN4dYzUgN1SBOBWkxwVPVnInI1cAYgwAJVXRL0yMJArbeOx1ftZcaI3pw6wlpvzAmGAJk+y1nArEb2m+1OGnsI+Kmqbm3Fa4N2lmtMW6nqBaGOoSm5pVUs3nCQr84Yan1vOjm/5qJS1ZeBl4McS9h5a8thsooq+dXlk0MdiglPjXXIangWugEYoaplInIp8Aow1s/XOiuDdJZrTDTyfLyPmro663tjmu5LIyKr3J+lIlLi8ygVkZKOCzE0VJV/rdjD6H5dOX9C/1CHY8JTFjDMZ3koTivNMapaoqpl7vO3gHgRSfHntcZEIhG5UkSygNnAm24Xhw5RUV3LM2v2M2/yQFJTbNybzq7JFhxVnev+7N5x4YSPNRmFbD1Uwv9ddRIxNneJadxaYKyIjAQOAtcBX/PdQUQGAjmqqiIyE+ekogA40tJrjYlEbheGkHRjeDXtECVVtdw2d2Qo3t6EGX/GwXnan3XR5uk1++iVHM+XTxkS6lBMmFLVWuC7wLvAdmChqm4VkTtE5A53t68AW9w+OI8A16mj0dd2/KcwJjqoKk+v3s/EQT04dUTvUIdjwoA/fXCO64AiInHAqcEJJzzklFTx7tYcbps7kqR4m3PKNM297PRWg3WP+Tz/O/B3f19rjGmbDQeOsC27hD9ceZKNV2aA5vvg3CUipcBU3/43QA7g922Akei5Tw9Qp8oNs+yOFWOMiQTPrNlP98Q45k8bHOpQTJhoMsFR1T+6/W8eUNUe7qO7qvZV1bva86Yicp+IHPQZ3fXSJvbr8JFea7x1PP/ZAc4e188mZzPGmAhQWF7Nm5uyufrUoXRN9OvmYNMJ+DMOzl0i0hvn1tYkn/XtnWfkIVX9c1MbQzXS63vbc8gtPcofZo0I5tsYY4wJkFc+P0i1t47rZ1qru/lCiwmOiHwT+AHObaxpwOnAauC8oEYWopFeX1qXxYAeiZxrt4YbY0xEeHlDFlOH9mT8wE55069pgj9zSv0AOA3Yr6rnAqcAeQF47++KyCYRecJtIWqosZFem7ylSURuF5F1IrIuL69t4eWWVvHhrjyumj6UWLs13Bhjwt62QyVsPVTC1dOHhjoUE2b8SXCqVLUKQEQSVXUHML6lF4nIchHZ0shjPvBPYDQwDcgG/tLYIRpZ1+Qorqq6QFVnqOqMfv36+fGxTvTK5wfx1ilfOdW+KMYYEwle3pBFfKxwxcnWudgcz5/eWFki0gtniPllIlKEHyOu+jtXiYj8G3ijsfelA0d6VVVeWpfF9OG9GN2vW7DexhhjTIDUeOt45fODXDBxgM0abk7gTyfjK92n94nIB0BP4J32vKmIDFLVbHfxSmBLI7u1OEpsIG3KKmZ3bhl/vOqkYL2FMcaYAPpodx4F5dV2eco0qskER0R6qGqJiPhOo73Z/dkNKGzH+94vItNwLjntA77lvudg4HFVvVRVa0WkfqTXWOCJYI70uuTzgyTExfClqYOC9RbGGGMC6I2N2fTsEs9Z49rWLcFEt+ZacJ4DLgPW4yQi0uBnm6dqVdWbmlh/CLjUZ7lDRnr11ilvbc7m3PH96JEUH+y3M8YY005VNV6Wbsvh0pMGkhDnT3dS09k0N9nmZe7PqJ+17LO9heSWHuVy66RmjDERYeWuPMqO1nLZVKu3TeP8mWzzVRG5XkSSOyKgUHhj0yG6xMdyno19Y4wxEeGNTdn0To5n9ui+oQ7FhCl/2vUeBM4EtovISyLyFRFJaulFkaLWW8fbWw5zwaQBJCfYEN/GGBPuqmq8LN+ew7wpg4iPtctTpnH+3EW1AljhTp1wHvA/wBNAjyDH1iE+2VNAYXk1l1nnYmOMiQgf7sylotpr9bZpll+pr4h0Aa4G7sAZ1fi/wQyqI721OZtuiXGcbb3wTRu0NCGsiNzgjti9SUQ+EZGTfbbtE5HN7oSz6zo2cmMi19KtOfRKjmfWyD4t72w6LX/monoRmIUz9s2jwIeqWhfswDqCt05Zvj2Hc8b3Iyk+NtThmAjj54Swe4GzVbVIRC4BFuB8n+qdq6r5HRa0MRGu1lvH+ztzOW98f+Ls8pRphj+dTp4Evqaq3mAH09HSMovIL6vmwkkDQh2KiUwtTgirqp/47L8GZ0RuY0wbrd9fxJGKGi6wetu0wJ/0dyVwl4gsABCRsSJyWXDD6hhLt+UQFyOcM97unjJt0qoJYYHbgLd9lhVYKiLrReT2pl4UiIlkjYkWy7fnkBAbY4P7mRb5k+A8CVQDc9zlLOB3QYuoAy3bmsPs0X3p2cUG9zNt4veEsCJyLk6C8wuf1Weo6nTgEuA7InJWY68NxESyxkQDVWXZthxOH92Xbol216tpnj8JzmhVvR+oAVDVShqv2CNKem4ZGfnldnnKtIdfE8KKyFTgcWC+qhbUr3dH7kZVc4ElOJe8jDFN2JNXzr6CCi6caK3upmX+JDjV7l1UCiAio4GjQY2qAyzddhiACyZagmPa7NiEsCKSgDMh7Gu+O4jIcGAxcJOq7vJZ31VEutc/By6i8UlnjTGu5dtzADjf6m3jB3/a+H6FcwfVMBF5FjgDuCWYQXWED3bkMnlwDwb36hLqUEyEampCWBG5w93+GHAv0Bf4h4gA1KrqDGAAsMRdFwc8p6rvhOBjGBMxVuzMY8LA7lZvG7/4M9DfMhHZAJyOc2nqB5F+W2txZQ0bDhzh22ePDnUoJsI1NiGsm9jUP/8m8M1GXpcBnNxwvTGmcRXVtazbX8itZ0T99IgmQJpMcERkeoNV2e7P4SIyXFU3BC+s4PokPR9vnXL2eOuwaYwxkWBNRgE1XuWssVZvG/8014LzF/dnEjAD2IjTgjMV+BSYG9zQgufDnXl0T4rjlGG9Qh2KMcYYP6zclU9SfAwzUnuHOhQTIZrsZKyq56rqucB+YLp7m+qpwClAekcFGGiqyopdeZw5NsVGwTTGmAixcnces0b2tVHnjd/8+Q8/QVU31y+o6hZgWtAiCrJdOWUcLqmyuaeMMSZCZBVVkJFXzpljU0Idiokg/txFtV1EHgeewblV/EZge1CjCqIVu3IBbBRMY4yJEB/tdu5rsRNT0xr+JDi3At8GfuAurwT+2Z43dSfwHO8u9gKOqOq0RvbbB5QCXr64vbZdVu7KZ/yA7gzqabcZGmNMJFiVns/AHkmM6d8t1KGYCOLPbeJVwEPuIyBU9dr65yLyF6C4md0DNtvy0Vova/cVcuPpIwJxOGOMMUGmqnyaUcCZY/vhjhtljF9COpmHOKX1GuC8jni/zw8c4WhtHbNH9e2ItzPGGNNOe/LKyC+r5vRRfUIdiokwob6N6EwgR1V3N7Hdr9mWwb8Zlz/ZU0CMwEz7ohhjTFCIyAMiskNENonIEhHp1Z7jrc4oBOB0OzE1rdRigiMiX/VnXSP7LBeRLY085vvsdj3wfDOH8Wu2ZfBvxuXVe/I5aUhPeiTZ7OHGGBMky4ApqjoV2AXc1Z6DrckoYFDPJIb3SQ5IcKbz8KcFp7HC2WKBVdULVHVKI49XAUQkDrgKeLGZYwRstuWK6lrSMo8we7TdZmiMMcGiqktVtdZdXAMMbcex+DSjkNNH9bX+N6bVmpuq4RLgUmCIiDzis6kHUNv4q1rlAmCHqmY18f5dgRhVLfWZbfk3bX2zdfuKqPEqc0ZbM6cxxnSQb9DMSazb9eB2gOHDh5+wfU9eOfllR63/jWmT5joZHwLWAVcA633WlwI/CsB7X0eDy1MiMhh4XFUvJcCzLX+yp4D4WLFhvo0xpp1EZDkwsJFNd/u00t+NczL8bFPHUdUFwAKAGTNmaMPtazIKAJg10k5MTes1meCo6kZgo4gsAcpV1QsgIrFAYnvfWFVvaWTdIZxWo4DPtrx6Tz7ThvUiOSGkN44ZY0zEU9ULmtsuIjcDlwHnq+oJiYu/1mQUMLBHEiP6Wv8b03r+/LdfinM5qcxd7uKumxOsoAJNVTkttQ+pKV1DHYoxxkQ1EZkH/AI4W1Ur2nOsyYN7MqZ/N+t/Y9rEnwQnSVXrkxtUtUxEIiqdFhHuuWxSqMMwUcitzP8KxOJcXv2/BtvF3X4pUAHcoqob/HmtMRHq7zit/MvcxGSNqt7RlgN9+5zRgYzLdDL+JDjlIjLdp1I+FagMbljGhD/3cu2jwIVAFrBWRF5T1W0+u10CjHUfs3CmOZnl52uNiTiqOibUMRgD/iU4PwReEpFD7vIg4Nqmdzem05gJpLv9xRCRF4D5gG+SMh94yu2HsEZEeonIICDVj9eeICOvnGv/tTrgH8QYY6KNP3NRrRWRCTiTYwrOrd01QY+sHdavX58vIvsb2ZQCBGReqyCzOAOvsVjbOynZECDTZzkLp5WmpX2G+Pla4PhbaYGjC++Ys6UdMQdKOP3tLZbGjW95l/Bh9XaHiZQ4oelY/aq7W0xwRCQJ+F9gLs7UCR+JyGPuJJxhSVUbHcpYRNYFYkbyYLM4Ay9IsTbW87HhHSNN7ePPa52VPrfShsvvPFziAIulKSKyLtQxtIbV2x0jUuKE9sfqzyWqp3DGvvmbu3w98DTQ4nQNxkS5LGCYz/JQnPGj/NknwY/XGmOMaSN/Epzxquo7Hs0HIrIxWAEZE0HWAmNFZCRwEGfwyq812Oc14LtuH5tZQLGqZotInh+vNcYY00b+JDifi8jpqroGQERmAR8HN6ygWRDqAPxkcQZewGNV1VoR+S7wLs6t3k+o6lYRucPd/hjwFs4t4uk4t4nf2txrQ/E52ihc4gCLpSnhFEt7RMrnsDgDr12xSkuDTIrIdpzOagfcVcOB7UAdoO6MscYYY4wxYcOfBKfZ3sqq2livd2OMMcaYkGkxwTHGGGOMiTQxoQ6gI4jIPBHZKSLpInJnGMTzhIjkisgWn3V9RGSZiOx2f/b22XaXG/tOEbm4A+McJiIfiMh2EdkqIj8Ix1hFJElEPhORjW6cvw7HOJvTUhkVxyPu9k0iMt3f1wYhlhvcGDaJyCcicrLPtn0isllE0gJxm7IfsZwjIsXu+6WJyL3+vjbAcfzMJ4YtIuIVkT7utoD9ThqrOxps77ByEmzhFq/V2wGPM/j1tqpG9QOnA+ceYBTOrbkbgUkhjuksYDqwxWfd/cCd7vM7gT+5zye5MScCI93PEttBcQ4CprvPuwO73HjCKlacMWW6uc/jgU+B08MtzvaUUZyOym+7n/V04NNglG8/Y5kD9HafX1Ifi7u8D0jpwN/LOcAbbXltIONosP/lwPtB+p2cUHeEopwE+xGO8Tb2uw/HOgart489OkMLzrHh9FW1GqgfEj9kVHUlUNhg9Xzgv+7z/wJf9ln/gqoeVdW9OHfjzOygOLPVnYNMVUtxOpcPCbdY1VE/IWy8+9Bwi7MZ/pTRY1M+qHNHY/2UD4Eu3y0eT1U/UdUid3ENzhg+wdCezxbI30trj3U98Hwb36tZTdQdvjqqnARb2MVr9XbA4wx6vd0ZEpymhsoPNwNUNRucAgr0d9eHRfwikgqcgpNlh12sIhIrImlALrBMVcMyzib4E09rpnxoz2dp7fFuw2kxqKfAUhFZL84UE+3hbyyz3Wbut0VkcitfG8g4EJFkYB7wss/qQP5OWtJR5STYIiXesK5jOnu97c84OJHO7yHxw1TI4xeRbjgV9g9VtUSksZCcXRtZ1yGxqqoXmCYivYAlIjKlmd1D/jttoEOmfAhgLM6OIufiJDhzfVafoaqHRKQ/sExEdrhnvsGKZQMwQlXLRORS4BWcmdsD+XtpzbEuBz5WVd8z/UD+TlrSUeUk2CIt3oZCHr/V252jBcef4fTDQY7blIz7M9ddH9L4RSQe50vyrKouDudYAVT1CPAhzll02MbZQHumfAj0Z/HreCIyFXgcmK+qBfXrVfWQ+zMXWEL7mrpbjEVVS+qbuVX1LSBeRFL8/RyBisPHdTS4PBXg30lLOqqcBFukxBuWdYzV218cOKofOK1UGTidkuo7q00Og7hSOb6z2gMc37Hqfvf5ZI7vWJVBx3VWE5y5yB5usD6sYgX6Ab3c512Aj4DLwi3O9pRR4Esc33n0s2CUbz9jGY5z/XtOg/Vdge4+zz8B5gU5loF8MdzFTJwBSSWQvxd/jwX0xOmj0TVYvxP3OKk03cm4Q8pJOHwnQhTXcb/7cKxjsHr7i/cIdYHpoD/4pTg9yfcAd4dBPM8D2UANTlZ6G9AXeA/Y7f7s47P/3W7sO4FLOjDO+hnkNwFp7uPScIsVmAp87sa5BbjXXR9Wcba2jAJ3AHe4zwV41N2+GZgRrPLtRyyPA0U+ZWKdu36UWwFtBLZ2UCzfdd9rI06H5znNvTZYcbjLt+B0gvR9XUB/JzRed4SknITiOxHieKzeDmycQa+3baA/Y4wxxkSdztAHxxhjjDGdjCU4xhhjjIk6luAYY4wxJupYgmOMMcaYqGMJjjHGGGOijiU4xhhjjIk6luAYY4wxJupYgmOMMcaYqGMJjjHGGGOijiU4xhhjjIk6luAYY4wxJupYgmOMMcaYqGMJTgQTERWRMaGOw5hoJCK/FJHHQx2HMaZtLMEJAhHZJyKVIlLm8/h7qOMKNRG5T0SeCXUcpu1E5Gsiss4t09ki8raIzA11XO0lIueISJbvOlX9g6p+M1QxmcgnIneJyFsN1u1uYt11HRtd9LMEJ3guV9VuPo/vhjogY9pDRH4MPAz8ARgADAf+AcwPYVjGhLOVwBkiEgsgIgOBeGB6g3Vj3H1NAFmC04FE5BYR+VhEHhKRIyKSISJz3PWZIpIrIjf77O8RkcdEZJmIlIrIChEZ0cSxe4rIUyKSJyL7ReQeEYkRkUQRKRSRk3z27e+2MPWrP3MVkZ+7758tIl8WkUtFZJf72l/6vDZGRO4UkT0iUiAiC0Wkj7st1b1sdrOIHBCRfBG52902D/glcK179r8xWL9nE3gi0hP4DfAdVV2squWqWqOqr6vqz9xy9rCIHHIfD4tIovva+jL2E58ydqvPsS8VkW1uGT8oIj91198iIqsaxHHssqz7/fiH24pU5n63BrrvXSQiO0TkFJ/X7nPPqLe5258UkSQR6Qq8DQz2aXEd3LDFUUSuEJGt7nf3QxGZ2ODYPxWRTSJSLCIvikhScP4aJoKsxUloprnLZwEfADsbrNsDXCwi293vQYaIfMv3QG4dne1+v77Z4LuQKCJ/duvdHPf/RpcO+HxhzRKcjjcL2AT0BZ4DXgBOw8ngbwT+LiLdfPa/AfgtkAKkAc82cdy/AT2BUcDZwNeBW1X1qPseN/rsez2wXFXz3OWBQBIwBLgX+Le7/6nAmcC9IjLK3ff7wJfd9xgMFAGPNohlLjAeON997URVfQfnzP9Ft0Xr5OZ+SSbszMYpI0ua2H43cDpOpX0yMBO4x2f7QJzyOQS4DXhURHq72/4DfEtVuwNTgPdbEdc17vukAEeB1cAGd3kR8GCD/W8ALgZGA+OAe1S1HLgEOOTT4nrI90UiMg54Hvgh0A94C3hdRBIaxDIPGAlMBW5pxecwUUhVq4FPcZIY3J8fAasarFsJ5AKXAT2AW4GHRGQ6HDtB/DFwAc7/irMbvNWfcMrzNHd7fV3euamqPQL8APYBZcARn8f/4FR4u332OwlQYIDPugJgmvvcA7zgs60b4AWGucuKU5hjcSr3ST77fgv40H0+C8gEYtzldcA17vNzgEog1l3u7h53ls+x1gNfdp9vB8732TYIqAHigFT3tUN9tn8GXOc+vw94JtR/H3u0qUzfABxuZvse4FKf5YuBfQ3KWJzP9lzgdPf5Abe89mhwzFuAVQ3WKTDGfe4B/u2z7XvAdp/lk4AjPsv7gDt8li8F9vjEmNXgvY6VV+D/AQt9tsUAB4FzfI59o8/2+4HHQv13s0foH245WuI+3wiMxUmEfdfd3MjrXgF+4D5/Avijz7YxfFH/C1AOjPbZPhvYG+rPHuqHteAEz5dVtZfP49/u+hyffSoBVLXhOt8WnMz6J6paBhTitJz4SgESgP0+6/bjZPGo6qc4X4CzRWQCzpfiNZ99C1TV6xtTI3HWxzQCWOI20x/BSXi8OH0y6h32eV7R4POYyFQApIhIXBPbB3Ni+fMtpwWqWuuz7FsursZJNva7l2FntyKuhuW0ue8S+HyfGomxOcd9PlWtc481xGcfK/emMSuBuW6LZT9V3Q18Asxx100BVorIJSKyxu0WcATnO5HiHmMwx5dd3+f9gGRgvU+9/I67vlOzBCf8Dat/4l666gMcarBPPk4rim//nOE4Z5j1/otz2ekmYJGqVrUxnkzgkgbJW5KqHmzxlc4Zh4lMq4EqnMuTjTnEieWvYTltlKquVdX5QH+cs9aF7qZynIobONYZs72G+Tz3jbGlsnnc5xMRcY/lT7k3ndtqnMuztwMfA6hqCU6Zut39eQh4GfgzTot+L5zLoOIeIxsY6nNM33Kcj5PMT/apk3uqaqdPsC3BCX+Xishc91r/b4FPVdU3e8dtfVkI/F5EuovTEfnHgO8t2U8DV+IkOU+1I57H3PcZASBOR2V/76LJAVJFxMpdhFHVYpxr+o+6ndCTRSTePeu8H6d/yj1ueUhx921xSAARSRCRG0Skp6rWACU4LYLgNN1PFpFpbofd+wLwUb4jIkPF6Rj/S+BFd30O0NftTN2YhcCXROR8EYkHfoJzWfiTAMRkopiqVuJ0C/gxTv+beqvcdStxWuATgTygVkQuAS7y2XchcKuITBSRZHz617itif/G6bPTH0BEhojIxcH7VJHB/tEEz+ty/Dg4TXXObMlzwK9wLk2ditMXojHfwznjzcD54jyHc90WAFXNwul8qRz/JWutv+Jc3loqIqXAGpw+Pv54yf1ZICIb2hGDCQFVfRCnQr4HpyLOBL6L0+ryO5xKfBOwGaes/c7PQ98E7BOREuAO3A7xqroL586t5cBunHLdXs8BS3G+Jxn1MarqDpwkLcNt5j/u0pWq7nTj+hvOGfPlOENBVAcgJhP9VuC0UPqW4Y/cdStVtRTnBo6FODdufA2fbgSq+jbwCM4dWOk4rULgJNkAv3DXr3G/R8txbvTo1MTtkGTCkIh4cDo+3tPSvn4e7wmcO0UCcjxjIomI7AO+qarLQx2LMe3hDlGwBUhs0LfN+LAWnE5CRFKBq3BuyTXGGBNBRORK95Jub5zbwl+35KZ5luB0AiLyW5xs/wFV3RvqeIwxxrTat3AuDe/B6af27dCGE/7sEpUxxhhjoo614BhjjDEm6jQ1aFdES0lJ0dTU1FCHYcLc+vXr81U1ogbDsrJt/BVp5dvKtvGXv2U7IhIcdx6Ov+JMSfC4qv5fc/unpqaybt26DonNRC4R2d/yXkGPwcq2CYpQl28r2yZY/C3bYX+JSpwp5R/FmQxvEnC9iEwKbVTGtJ+VbROtrGybcBD2CQ7OrMTpqprhDqr1AtDsyLkFBQWkpaUB4PV68Xg8bNq0CYCamho8Hg9btmwBoKqqCo/Hw/bt2wGoqKjA4/Gwc+dOAMrKyvB4PKSnpwNQXFyMx+MhIyMDgKKiIjweD/v27QMgPz8fj8dDZqYz2HBubi4ej4eDB50R3Q8fPozH4+HwYWfamoMHD+LxeMjNzQUgMzMTj8dDfn4+APv27cPj8VBUVARARkYGHo+H4uJiANLT0/F4PJSVlQGwc+dOPB4PFRUVAGzfvh2Px0NVlTMzw5YtW/B4PNTU1ACwadMmPB4PXq8zeGxaWhoej+fY73L9+vU89dQXAx+vXbuWZ5/9YkLzNWvW8Pzzzx9b/uSTT1i4cOGx5VWrVrFo0aJjyx9+uIKFi17mSEU1uaVVvPb2Up55cRG7c0rZeqiY55a8yRPPLeKTPfms2JXHv59/hX8+/RKvbTzEq2kHefTpRfz9qUUsWp/FwnWZ/PXJF3nkqUU8++l+nlmznwcff46H//syT368l/+s2ssDC57hb8+8cuz9Fy9ezIoVK04oMyFiZTuKyvaKFSt4qUHZfvbFRaTnlrLtUAnPL3mTJ59fxOo9BXy0O4/HX3DK9hubDvHK5wf5+1NO2X5pXSYL1zZdtj0f7+UJt2w//NQX44da2bayXS9Y9XZR+fH1dn3Zrq+3V+8pYOWuPB536+3XNzplu77efmldZrP1tm/Zfui/S6jx1gFtL9uRcIlqCMdPLJZFIyPnisjtOPN6MGTIkIabTTtV1XgpKK8mu7iKoopq3tlymJKqGg5lFFCWX8avXt1CSVUtkptJbGURLz32CZU1XvqW7aVLbSm//81SKqq9TCCLnlLFz9ctA+CUuIN0lWru+XwlAKfGZZEotfxm06cAnBZ3mFipY83WzwGYGZ8HwGfbNgJwenwBXo1h7Tan4psTX8RRjWP99m0AzI0/gsYf5Xsd96tqDSvbYaDWW0dhRTU5JVUUV1azfFsOJVU1ZO0tpDSvnN++sY2SyhrqcrKIrShk0b9WU1XjpW9pBom1pfzxt8uoqPYyTjPpKVX8rEHZvrtB2f71xjVAc2Xb+afuT9mulEp+2GG/qVaxsh0Gvqi3Kzuo3g5c2S7XCr7lrSM+tu3tMGF/m7iIfBW4WFW/6S7fBMxU1Sb/Z82YMUPtWq5/VJWiihqyiirILKwkq6iCnJKj5JUdJa+0irzSo+SVHqWkqunxpESge2Ic3ZPi6ZYYR1JCLF3iY0hOiKNLfCxdEmLpEh9LckIsSe5yYlwM8bExJMTGEBcrxMe6y3FfPHcecuyniBAjQoxAjAji/qx/7rscIyDH5qmDnsnxjcQt61V1RlB+sX6wsh185UdrySqqJLOwgqyiCg4VV5FfWl++j5JfdpSC8mqaqwa7JsTSo4tTtuvLcHKCU46PPY93HwlxJMbFkBDnW3aPfx4XK265d9YnxMYQGyNNlmXfZeGL9QCK0is5odG4Q1m+rWwHl6pSWF5NZpFTZ2cVVXK4uIq8sqPHle/SVtTb9fV0l4Qv6uym6u3jy/bxywlxQlzMF+vbU7Z7dolHRBqJ3b+yHQktOFkcP3PqUPycpdh8obq2joz8MnbllLE7p5RdOaXsL6ggs7CC8mrvcft2TYilX/dE+nVPZPzA7swdk0K/7omkdEukV3I8PZLi6dElnp5dnJ/dE+OIiTmxEJoWWdkOgLo6JbOogt05ZezKLSU9p4w9eWVkFlVSWH78VFEJcTH06+aU7aG9kzlleG+nrHdLoHfXBHrWl2u3jPdIiiOuHWeQnZiV7QCoqvGSkVfO7lynzt6dU8a+gnKyiiqpaFBvd0uMo79bT08c2IMzxyT41NsJ9OgSR4+kzlVvR0KCsxYYKyIjgYPAdTgTkZkm1Hrr2JlTSlrmEdIOHGFj1hH25JXjrXNOU2NjhBF9kxmV0pXTR/VlaO8uDOuTzLDeyQzt04UeSSe2dpigsLLdSqpKZmEln2cWsTGzmLTMIrZll1BVU3dsn4E9khjTvxsXT+7JsD5dGNo7mWG9nZ8p3RIaPSM0AWdlu5VqvHXsyC4lLeuLejsjrwy32iY2RhiZ0pWRKV2ZO6bfsXp7aO8uDO3dhe5Wb58g7BMcVa0Vke8C7+LcbviEqm4NcVhhpa5O2ZZdwqr0fFbtzmf9/iIqa5zsvndyPNOG9eKiSQMZO6Ab4wZ0Z1S/riTGxYY4amNl2z+ZhRV8nJ7PqvR8Vu8poMBtlUmMi+GkIT25fuZwxg/oztgB3Rk7oJsl6GHAynbLvHXK5oPFfJyez8fpTr19tNZJ1Pt2TWDasF5cetIgxg3oxtj+3RmZ0pWEOGtNbI2wT3AAVPUt4K1QxxFOqmq8rNiVx9ubs1m5O/9YU/y4Ad249rRhnDK8F9OG9WJ4n2Q7Yw1jVrZPVFenbDhQxNtbDrN8ew77C5w7Swb0SOTscf2YkdqHk4f1ZNyA7u3qgGiCy8r2icqP1vLBzlze2XKYlbvyjvVtnDioBzfMGsH0EU69PaRXF6u3AyAiEhzjqPHW8cGOXF7deIgPduRSUe2ld3I8547vz9yxKcwdk0L/HkmhDtOYVlNVNmUV8/KGLN7ecpi80qMkxMYwd2wKt85JZe7YFEb362aVvok41bV1LN+ew5LPD7JyVx5Ha+vo2zWBeVMGMndsP+aM7ktKt8RQhxmVLMGJAPsLynlhbSaL1meRV3qUlG4JXHnKEC49aRCzRvaxTpAmYhVX1PDyBmdMox2HS0mMi+G8Cf2ZN2Ug503ob/0KTMRKzy3jhc8OsPjzgxSWVzOgRyLXzxzOvCkDOS21D7FR3sE3HFiCE8bSMo/w2Id7eHfbYQQ4b0J/rj1tOOeO72dJjYloB49U8p+P9vLC2gNUVHuZOrQnv/vyFK6YNtj60JiIpaqs21/Ev1bsYfn2XOJihAsnDeCa04Zx1th+ltR0MEtwwtBnewt5cNlO1mQU0iMpjv89ZzQ3nZ7KwJ52+clEtszCCh5avotX0w4hwOUnD+a2uSOZMqRnqEMzpl0+Ts/nL0t3suHAEXonx/OD88dy4+kj6NfdLj+FiiU4YWRXTil/ensH7+3IpX/3RO6+dCLXzxpOt0T7M5nIVlB2lL+9n86zn+4nRoRb5qRy29yRDO7VJdShGdMuWw4W86d3dvDR7nwG90ziN/Mn89VTh9Elwe5UDTX7zxkGyo7W8pelO/nvJ/vomhjHz+eN59Y5I+0LYiJeXZ3y7GcHuP+dHVRUe7lmxjB+cP5Ya400Ea+4soYH3t3Bs58eoGeXeO750kRuPH0ESfFWb4cLS3BCbOnWw/zqta0cLqnixlkj+PGF4+jdtfGh142JJNuzS7hr8WbSMo8wZ3RffjN/CmP6dwt1WMa0i6ry5uZsfv36NgrKjnLLnFR+dOE46zsWhizBCZHyo7Xc++pWXt6QxYSB3Xn0hulMH9471GEZ0251dcoTH+/lT+/soEdSPA9dezJfnjbEbvE2Ea+kqoZ7lmzhtY2HOGlIT5685TTrPxbGLMEJgS0Hi/ne85+zv6Cc7583hu+dP9YGLDNRIbe0ip8s3MhHu/O5aNIA/nT1VGuRNFFhw4Eivv/852QXV/HTi8bx7XPG2F1RYc4SnA62cG0md7+ymb5dE3nuf07n9FF9Qx2SMQHx+YEibn96PaVVNfz+yil8beZwa7UxUeGZNfu577WtDOyZxMJvzebUEdbaHgkswekg3jrlj29t5/FVezlzbAqPXHeKndmaqLF4QxZ3Lt7MgB6JPH3bGUwY2CPUIRnTbjXeOn77xjaeWr2fc8f346/Xn2J9bSKIJTgdoKrGy3ef+5zl23O4ZU4q93xpog3UZ6KCqvK399N5cNkuTh/Vh3/ecKol7iYqVFZ7ueOZ9azYlce3zhrFz+dNsEtSEcYSnCArO1rL//x3HWv2FvCb+ZP5+uzUUIdkTEDU1Sm/e3M7T3y8l6umD+FPV0+1vmQmKpRU1XCbZy3r9xfxx6tO4vqZw0MdkmkDS3CCqLiyhpuf+IzNB4t56JppfPmUIaEOyZiAqKtT7ly8iYXrsrhlTir3XjaJGDu7NVHgSEU1N/7nU3YeLuVv10/nS1MHhTok00aW4ARJRXUttz75GVsPFfPPG6Zz0eSBoQ7JmIBQVX712lYWrsvi++eN4UcXjrPOxCYqlB2t5eYn17LrcBkLvj6Dc8f3D3VIph0swQmCo7VevvX0etIyj/APS25MFFFV/vTOTp5es5/bzxplyY2JGlU1Xm5/ah1bDjonpZbcRD5LcAKsrk75sTsOyP1XT2XeFGveNNHjP6v28tiKPdwwazh3XTLBkhsTFerqlB+9mMYnewp46NqT7aQ0SliPwAD763u7eXNTNndeMoFrThsW6nCMCZj3d+Twh7e2M2/yQH47f4olNyZqPLR8F29vOcw9X5rIlacMDXU4JkAswQmgNzdl89f3dvOVU4fyrbNGhTocYwJm5+FSvv98GpMG9+DBa0+2DsUmaryadpC/vZ/OtTOGcdvckaEOxwSQJTgBsu1QCT95KY1TR/Tm91fa2a1pnog8ICI7RGSTiCwRkV6hjqkpxZU1fPOptSQnxPL4108jOcGubJvosPVQMT9ftImZqX347Zet3o42luAEQEV1Ld99fgM9kuJ57MZTSYyLDXVIJvwtA6ao6lRgF3BXiONplKpy1+JNZB+p4p83nsrAnkmhDsmYgCg/Wsv3nvucXsnx/PPG6STE2b/DaGN/0QC477Wt7M0v5+HrptGve2KowzERQFWXqmqtu7gGCMsL/898eoC3Nh/mpxePt/l3TFS599Wt7C0o5+FrT6FvN6u3o5ElOO302sZDLFyXxXfOGcOc0SmhDsdEpm8Abze1UURuF5F1IrIuLy+vw4Lanl3Cb9/Yxtnj+nH7mdanzESPVz4/yMsbsvjeuWOYPdomPI5WdjG9HXJKqrh7yWamD+/FDy4YG+pwTJgRkeVAY/eb3q2qr7r73A3UAs82dRxVXQAsAJgxY4YGIdQT1Hjr+OlLG+mRFM+D11inYhM9ckqq+H+vbmHGiN58/3yrt6OZJTjtcO+rW6iureMv10yzOXjMCVT1gua2i8jNwGXA+araIYmLvxaszGDroRIeu3G6Nd+bqFJfbz/w1ZNt0uMoZ3/dNnp7czbvbs3hhxeMY2RK11CHYyKMiMwDfgFcoaoVoY7HV3puKX9dvpsvnTTIBqo0UcXq7c7FEpw2KK6o4d7XtjJ5cA/+50wbN8G0yd+B7sAyEUkTkcdCHRA4I7r+fNEmkhNjue+KyaEOx5iAqa+3pwyxeruzsEtUbfCXZTspKDvKk7ecZk2cpk1UdUyoY2jM4s8PsuHAEf781ZPtjkATVR5+b5fV252M/ZVbaXdOKc9+eoAbZo1gypCeoQ7HmIApO1rLn97ZwbRhvbjqlCGhDseYgEnPLePp1fu5fuZwq7c7EUtwWul3b24nOSGWH104LtShGBNQj36QTl7pUX51+SS7a8pElT+8tZ0u8bH82OrtTsUSnFb4YGcuK3bl8YPzx9Kna0KowzEmYPYXlPOfj/Zy1fQhnDLcBvQz0WPFrjze35HL984fY3cEdjKW4Pip1lvH79/cTmrfZL4+OzXU4RgTUH9euou4WOEX8yaEOhRjAqauTvnDm9sZ0TeZm+ekhjoc08EswfHTaxsPkZ5bxi/mTbA5S0xU2XG4hDc2HeKWOakM6GFzTZn2EZGvishWEakTkRmhjOWNzdnszCnlJxeNtzkCO6Gw/k8dLjMu13rr+Ot7u5k0qAcXT25sYFpjItdDy3bRLSGO28+y6RhMQGwBrgJWhjIIb53y8PJdjB/QnctOsvGcOqOwTnAIkxmXF39+kP0FFfzownHW+dJElS0Hi3l3aw63nTmSXsnWr8y0n6puV9WdoY7jtY0Hycgr50cXjrV6u5MK6wQnHGZcrvHW8ch7u5k6tCcXTOzf0W9vTFA9uGwXPbvE8425NvCZiR7eOuWR99KZNKgHF02yVvfOKqwTnAZCMuPyG5sOkVVUyQ/OH4uInQWY6LHlYDHv78jl9rNG0SMpPtThmAgiIstFZEsjj/mtPE5Q6u13tx5mb3453ztvjLXedGIhH8k4nGdcVlX+tSKDsf27ce54a70x0WXBygy6JcZx4+kjQh2KiTAtTSTbiuMEqd7eQ2rfZC6yPpOdWsgTnHCecXnl7nx2HC7lga9MtbMAE1UyCyt4c3M23zgjlZ5drPXGRI9P9xayMauY3315CrFWb3dqYX2JKtQzLi9YuYcBPRKZP82GrTfR5T+r9iJgfW9MwInIlSKSBcwG3hSRdzvy/ReszKBv1wS+cmqHd9k0YSasExxCOOPy1kPFfJxewK1njLRxb0xUKa6o4cW1mVwxbTCDenYJdTgmyqjqElUdqqqJqjpAVS/uqPfek1fG+zty+frsVJLibdybzi7kl6iaE8oZl59Zs5+k+BiuP214qEIwJiheWp9JZY2Xb861cW9MdHlmzX7iY4WvzbJ624R/C05IFFfW8Mrnh5h/8hB6Jlv/BBM96uqUZ9bsZ8aI3kwa3CPU4RgTMBXVtSxan8UlUwbRr7vNOWUswWnU4g1ZVNZ4uWm23V1iosuq9Hz2FVRY2TZR57W0Q5RW1VrZNsdYgtOAqvL0mv2cMrwXU4b0DHU4xgTU02v207drAvOm2O2zJnqoKk+t3s+Egd2ZMaJ3qMMxYcISnAZWZxSQkVfOTTY2iIkyh45U8t72HK6bOcwmHjRRZVNWMduyS7jh9BE2IKs5xhKcBhaty6J7UhyX2uRsJsosWp9FncJ11nHeRJlF67NIjIth/rTBoQ7FhBFLcHyUVtXw1pZsLj95sN1iaDqEiPxURFREUoL5PqrKovVZzBndl2F9koP5VsZ0qKoaL69tPMS8KQNtyhFzHEtwfLy9+TBVNXU2QJTpECIyDLgQOBDs91q7r4gDhRVWtk3UeW97LsWVNVw93cq2OZ4lOD4Wrc9iVL+unDKsV6hDMZ3DQ8DPgaBPQfLSuky6JsRa52ITdV7ekMXAHkmcMSaojaAmAlmC49qXX85n+wr5yqlDrZOaCToRuQI4qKobg/1e5UdreXNzNl+aOojkhLAe29OYVsktrWLFrjyumj7E5p0yJ7DazvVK2kFE4MpTbN4pExgishxorMnkbuCXwEV+Hud24HaA4cNb30F46bbDVFR7rQnfRJ03NmbjrVOumm71tjmRJTg4HTDf2JTNzNQ+NjePCRhVvaCx9SJyEjAS2Oi2Fg4FNojITFU93MhxFgALAGbMmNHqy1mvb8xmcM8kTkvt09qXGhPW3tyczYSB3RnTv3uoQzFhyC5RATtzSknPLeOyk+0WQxN8qrpZVfuraqqqpgJZwPTGkpv2Kq6o4aPdeXxp6iBirAnfRJFDRypZv7+Iy63eNk2wBAenmTNGYN5k64Bposu7Ww9T41X7J2CizlubswFszDLTpE5/icq5PHWI2aP72gRtJiTcVpygeH3TIYb3SeYkm3bERJnXN2UzeXAPRqZ0DXUoJkx1+hacrYdK2FdQwWVT7QzXRJeCsqN8sqeAy6YOsjsDTVTJLKxgY+YRq7dNszp9gvPm5mxiY8QuT5mo887Ww3jr1P4JmKjzzhanu9qX7PKUaUanT3CWbcth1sg+9O6aEOpQjAmopVtzGNE3mYmD7A4TE12Wbc9hwsDuDO9r046YpnXqBGdvfjnpuWVcOGlAqEMxJqDKjtayek8BF04cYJenTFQpKq9m3b5Cq7dNizp1grN8Ww4AF0y0L4qJLit25lHtrbN/AibqfLAzlzq1etu0rFMnOMu2Oc2cNruyiTbLth2md3I8p47oHepQjAmo5dtz6N890e4MNC3qtAlOYXk16/YXcpGd4ZooU+Ot4/0duZw3YQBxsZ32K26i0NFaLyt25nH+xAE2cKVpUaet/d7f4TZzWoJjosxnewspqaq1y1Mm6qzJKKS82msnpsYvnTjByWFAD2vmNNFn2bYcEuJiOGtcSqhDMSag3t+eQ1J8DLNH9w11KCYCdMoEp9Zbx0e78zl7XD+7w8REnZW78pg9qi/JCZ1+oHITZVbuzmf2qL4kxceGOhQTATplgrMx6wilVbWcNa5fqEMxJqAyCyvIyC+3sm2iTmZhBXvzyzlzrJVt459OmeCs2JVPjMDcMdaEb6LLil15AJxtCY6JMit3O2Xbknfjr06Z4KzclcfUob3olWyjF5vosnJXHkN6dWF0P5uA0ESXlbvyGNwzycq28VunS3CKyqvZlHXEzgJM1Knx1vHJngLOsr5lJsrUeuv4JN3KtmmdTpfgrErPp07hbLvDxESZDfuLKDtaa5enTNRJyzxC6dFa639jWqXTJTgrd+XRPSmOk4f2CnUoxgTUil15xMUIc8bYLbQmuqzc7fSbPMPKtmmFTpXgqCqf7CngjNEpNsKriTofp+czbVgveiTFhzoUYwJq9Z58ThrS0/pNmlbpVP/ls4oqOXik0gaJMlGnpKqGzQeLmWNl20SZymovaZlHON3KtmmlTpXgrN5TAMDpo+yLYqLLZxmF1Cn2T8CEnIg8ICI7RGSTiCwRkV7tOd6GA0XUeNXqbdNqEZHgiMhPRURFpF09g9dkFNCnawLjBnQLVGjGhIXVGQUkxMUwfbjNHm5CbhkwRVWnAruAu9pzsE8zCogRmDHCyrZpnbBPcERkGHAhcKA9x1FV1mQUcPqoPnaboQkLIvI9EdkpIltF5P72HGv1ngJOHd7bhrA3IaeqS1W11l1cAwxtz/HWZBRy0pCedLe+ZaaVwj7BAR4Cfg5oew5yoLCCQ8VVzLZmThMGRORcYD4wVVUnA39u67GOVFSz/XCJ9S0z4egbwNtNbRSR20VknYisy8vLO2H7sf43Vm+bNgjrBEdErgAOqupGP/Zt9ouyJsP635iw8m3g/1T1KICq5rb1QGsyClHFEhzTYURkuYhsaeQx32efu4Fa4NmmjqOqC1R1hqrO6NfvxDFuPj9QRLW3zupt0yYhn25YRJYDAxvZdDfwS+Aif46jqguABQAzZsw4obVn9Z4CUrolMKa/9b8xYWEccKaI/B6oAn6qqmsb21FEbgduBxg+fPgJ29dkFNAlPtbGdjIdRlUvaG67iNwMXAacr6ptbn1fs7fQ6X+Tav1vTOuFPMFp6osiIicBI4GNbp+ZocAGEZmpqodb+R6syShk1qi+1v/GdJgWkvc4oDdwOnAasFBERjX2z8Cf5H1Gam8S4sK6QdZ0EiIyD/gFcLaqVrTnWGsyCqz/jWmzkCc4TVHVzUD/+mUR2QfMUNX8thzvqdtm0vbzCGNar7mzXBH5NrDYTWg+E5E6IAU48fpqC/5543Qqqr1tD9SYwPo7kAgsc08o16jqHW050J+/cjJFFdWBjM10ImGb4ASSiDBuQPdQh2GMr1eA84APRWQckAC0KXkf1c8uu5rwoapjAnWs4X2TGd43OVCHM51MxCQ4qpoa6hiMCaAngCdEZAtQDdzcnr4KxhhjjhcxCY4x0URVq4EbQx2HMcZEK4nGk0YRyQP2N7IphTZeBuhgFmfgNRbrCFU98d7UMGZlu8NESpzQdKwRVb6tbHeYSIkT2lm2ozLBaYqIrFPVGaGOoyUWZ+BFUqxtESmfz+IMvEiKtS0i5fNZnIHX3ljtvlJjjDHGRB1LcIwxxhgTdTpbgrMg1AH4yeIMvEiKtS0i5fNZnIEXSbG2RaR8Posz8NoVa6fqg2OMMcaYzqGzteAYY4wxphOwBMcYY4wxUadTJDgiMk9EdopIuojcGQbxPCEiue4otvXr+ojIMhHZ7f7s7bPtLjf2nSJycQfGOUxEPhCR7SKyVUR+EI6xikiSiHwmIhvdOH8djnEGg5XtNsdpZTvMWdluc5xWtuupalQ/gFhgDzAKZ76fjcCkEMd0FjAd2OKz7n7gTvf5ncCf3OeT3JgTcWZX3wPEdlCcg4Dp7vPuwC43nrCKFRCgm/s8HvgUZ5busIozCJ/bynbb47SyHcYPK9tWtgMRZ2dowZkJpKtqhjrD478AzA9lQKq6EihssHo+8F/3+X+BL/usf0FVj6rqXiAd5zN1RJzZqrrBfV4KbAeGhFus6ihzF+Pdh4ZbnEFgZbvtcVrZDm9Wttsep5VtV2dIcIYAmT7LWe66cDNAVbPBKaBAf3d9WMQvIqnAKThZdtjFKiKxIpIG5ALLVDUs4wywSPkcYf13sLIdliLlc4T136Gzl+3OkOBII+si6d74kMcvIt2Al4EfqmpJc7s2sq5DYlVVr6pOA4YCM0VkSjO7h/x3GiCR/jlCHr+V7bAV6Z8j5PFb2e4cCU4WMMxneShwKESxNCdHRAYBuD9z3fUhjV9E4nG+JM+q6uJwjhVAVY8AHwLzCOM4AyRSPkdY/h2sbIe1SPkcYfl3sLLt6AwJzlpgrIiMFJEE4DrgtRDH1JjXgJvd5zcDr/qsv05EEkVkJDAW+KwjAhIRAf4DbFfVB8M1VhHpJyK93OddgAuAHeEWZxBY2W4jK9thz8p2G1nZ9hHsntLh8AAuxelJvge4OwzieR7IBmpwstLbgL7Ae8Bu92cfn/3vdmPfCVzSgXHOxWkC3ASkuY9Lwy1WYCrwuRvnFuBed31YxRmkz25lu21xWtkO84eVbSvb7Y3TpmowxhhjTNTpDJeojDHGGNPJWIJjjDHGmKhjCY4xxhhjoo4lOMYYY4yJOpbgGGOMMSbqWIJjjDHGmKhjCY4xxhhjos7/B+GHCedkF3utAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 576x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 3, figsize=(8, 4))\n",
    "ax = axes.flatten()\n",
    "\n",
    "ax[0].plot(100*((irf['N'] + ss['N'])/ss['N']-1))\n",
    "ax[0].axhline(0, color='gray', linestyle=':')\n",
    "ax[0].set_title('Mass of establishments')\n",
    "ax[0].set_ylabel('pct deviation from ss')\n",
    "\n",
    "ax[1].axhline(0)\n",
    "ax[1].set_title('Markup')\n",
    "\n",
    "ax[2].plot(100*((tfp)/tfp[-1]-1))\n",
    "ax[2].axhline(0, color='gray', linestyle=':')\n",
    "ax[2].set_title('TFP')\n",
    "\n",
    "ax[3].plot(100*((irf['L'] + ss['L'])/ss['L']-1))\n",
    "ax[3].axhline(0, color='gray', linestyle=':')\n",
    "ax[3].set_title('Employment')\n",
    "\n",
    "ax[4].plot(100*((irf['C'] + ss['C'])/ss['C']-1))\n",
    "ax[4].axhline(0, color='gray', linestyle=':')\n",
    "ax[4].set_title('Consumption')\n",
    "\n",
    "ax[5].plot(100*((irf['w'] + ss['w'])/ss['w']-1))\n",
    "ax[5].axhline(0, color='gray', linestyle=':')\n",
    "ax[5].set_title('Wage')\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
